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ABSTRACT 

Supernovae (SN), the most energetic stellar feedback mechanism, are crucial for reg¬ 
ulating the interstellar medium (ISM) and launching galactic winds. We explore how 
supernova remnants (SNRs) create a multiphase medium by performing 3D hydrody- 
namical simulations at various SN rates, S, and ISM average densities, n. The evolution 
of a SNR in a self-consistently generated three-phase ISM is qualitatively different from 
that in a uniform or a two-phase warm/cold medium. By travelling faster and further 
in the low-density hot phase, the domain of a SNR increases by > 10 2 ' 5 . Varying n and 
S, we find that a steady state can only be achieved when the hot gas volume fraction 
,/V,hot ;$ 0.6 ± 0.1. Above that level, overlapping SNRs render connecting topology of 
the hot gas, and the ISM is subjected to thermal runaway. Photoelectric heating (PEH) 
has a surprisingly strong impact on /v,hot- For n > 3 cm -3 , a reasonable PEH rate is 
able to suppress the thermal runaway. Overall, we determine the critical SN rate for the 
onset of thermal runaway to be S cr it = 200(n/lcm 3 ) fc (E S N/10 51 erg) 1 kpc 3 Myr 3 , 
where k = (1.2, 2.7) for n < 1 and > 1 cm -3 , respectively. We present a fitting formula 
of the ISM pressure P(n, S), which can be used as an effective equation of state in 
cosmological simulations. Despite the 5 orders of magnitude span of (n, S), the average 
Mach number varies little: M ~ 0.5 ± 0.2, 1.2 ± 0.3, 2.3 ± 0.9 for the hot, warm and 
cold phases, respectively. 


1. Introduction 


The ISM permeates throughout galaxies, absorbs and scatters the light from astronomical 
objects, gives birth to stars, and provides material for galactic outflows. Knowledge about the 
ISM determines our ability to interpret astronomical data and to apprehend many fundamental 
astrophysical processes. The ISM is multiphase in nature. The hot component (T > 10 6 K) 


contains little mass, but occupies roughly half the volume in the solar neighbourhood (Snowden 
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et al. 1998), and may have the potential to generate galactic winds (e.g. Li & Wang 2013). How 
the multiphase structure forms and evolves, however, is still poorly understood. 

SN play a key role in shaping the ISM by injecting energy, momentum, mass and metals. They 
are the most energetic of the stellar feedback processes (radiation from stars yield more energy in 
total, but the efficiency of its coupling to gas is lower). SN explosions drive blast waves into the 
ISM, create the hot phase, generate turbulence and accelerate cosmic rays. A better understanding 
of SN feedback is therefore crucial to study the evolution of the multiphase ISM. 

Recent progress in observation and simulation makes it an ideal time to revisit the problem 
of SN feedback and the multiphase ISM. New observational techniques have started to resolve the 
individual star formation sites in external galaxies, allowing us to study stellar feedback in different 
environment at sub-kpc scale. A correlation between the gas energy and star formation 30 — 40 Myr 
ago, which is the typical lifetime of the progenitors of core collapse SN, implies the importance of 
SN feedback ( Stilp et al.|2013 ). Studies of 30 Doradus, a well-resolved star forming site in the Large 


Magellanic Cloud, offer precious information on stellar feedback in detail (e.g. Chu Sz Kennicutt 


1994, 

Pellegrini et al. 

2010 

Lopez et al. 

2011 

). Light echoes from SN themselves reveal the 3D 

structure of the surrounding ISM ( 

Xu et al. 

1995 

Kim et al. 

2008 

). The velocity dispersion of 


Ha emission may have a mild correlation with the local star formation rate at both low and high 
redshifts (Green et al. 2010 2014 Genzel et al.||20fl Arribas et al. 2014 Moiseev et al. 2015). 


Ab initio cosmological simulations find that feedback is central in reproducing reasonable 
properties of galaxies. Simulations lacking feedback produce galaxies inconsistent with observations: 
the masses are too big, and baryons settle to the center and form big bulges, rather than producing 
extended thin disks (Navarro Sz Benz 1991 Navarro & White 1994| |Abadi et al. 2003); galactic 


winds, which are observed ubiquitously in the star-forming galaxies (e.g. Heckman et al. 2000; 


Steidel et al. 

2004 

2010 

Strickland & Heckman 

2009, 

Bouche et al. 

2012; 

Martin et al. 2012 

), are missing. 


By implementing “feedback” the above 
problems are largely mitigated, since feedback heats and ejects gas, preventing gas from turning 
into stars, launching galactic scale winds and fountains, which reduce the mass of the bulge and 


transport gas into outer parts of galaxies (£ 

jomrner-Larsen et al 

1999; Scannapieco et al. 2012 

Hummels Sz Bryan 

2012, 

Hopkins et al. 

2012 

Stinson et al. 2013; 

Hopkins et al. 2014 

Ubler et al. 

2014). However, due to the limited resolution, cosmological simulations can only implement SN 


feedback in a sub-grid manner, and most of the current simulations use ad hoc models which involve 
turning off standard physics such as cooling or gas dynamics. As a result, the predictive power of 
such simulations is seriously limited. A more physically-based model of SN feedback is essential to 
improve the current galaxy simulations, and to achieve this, a better understanding of the interplay 
between SN and the multiphase ISM is indispensable. 

The evolution of a SNR in a uniform medium is well-studied by many authors. Early ID 


simulations (e.g. Chevalier 1974. CiofR et al. 1988) showed that the evolution is characterized by 


four stages: free expansion, Sedov-Taylor, “pressure-driven snowplow” and momentum-conserving. 









































































































































- 3 - 


Thornton et al. (1998) revealed the feedback at later stages is dominated by kinetic energy while 


thermal feedback is negligible due to rapid cooling. Taking into account the kinetic feedback may 


have big implications for galaxy formation (Simpson et al. 2015, Nunez et al. 2015, in prep) 


SNRs evolve in a qualitatively different way when exploding in a clumpy medium. Pioneer 


work by Cowie et al. (1981) studied the remnant evolution and energy budget in a three phase 


medium with most volume hot, and found that the blast wave travels much faster and further 
in the low-density, hot medium, and the loss of explosion energy happens later compared to the 
case when the medium is uniform. Recent 3D simulations by Martizzi et al. (2015) and Kim &; 


Ostriker (2015) confirmed that the evolution and energy feedback of a SNR in a clumpy medium 


is quite different from the uniform case, but the terminal radial momentum is affected little by the 


inhomogeneity of the ISM. Walch &; Naab (2015) found similar results for molecular clouds, and 


the pre-SN ionization can slightly enhance the final momentum per SN. 


Recently works on SN feedback have mainly considered a relatively dense medium and neglect 
the third phase of the ISM - the tenuous hot gas. However, a medium with a significant fraction 
of volume hot exists widely, e.g. at the disk-halo interface and the galactic winds launching sites. 
Moreover, the hot phase being the most tenuous transports energy/momentum much faster (speed 
of sound/blast wave oc p ~ 0,5 ) and cools very inefficiently (cooling time oc /K 1 ). It is therefore very 
interesting to investigate the SN feedback in a three-phase ISM. 


The multiphase ISM adopted in the studies of SN evolution has been largely set up by hand. 
However, we should note that while SNRs are affected by the lumpiness of the ISM, they are also 
the major player in shaping the overall properties of the gas. Thus it is very important to study 
SN feedback self-consistently. McKee &; Ostriker (1977) proposed the underlying theory of the 


three-phase ISM regulated by SN. Their basic ideas were the following: 


(i) If the SN rate is sufficient to render the porosity of the ISM ~ 1, then the ISM is inevitably 
multi-phase with the majority of volume being hot; 

(ii) All three phases are in approximate pressure equilibrium (Spitzer 1956); 

(iii) Different phases are in dynamical equilibrium in terms of mass and energy: SN keep 
creating hot bubbles and compressing diffuse gas into dense clouds, while clouds are evaporated by 
the hot medium; energy input from SN balances dissipation in all phases. 


(Spitzer 1956 


These principles yield a complete set of equations, and the solutions determine the pressure 
and mass/volume quota for each phase, given the mean gas density and SN rate. Their results were 
in rough agreement with observations of the ISM in the solar neighbourhood. Despite some imper¬ 
fectness of the theory (for example, underestimating of the mass of warm gas, not considering the 
clustering of SN, etc), their main ideas are still widely accepted as the standard picture of the local 
ISM. Numerical simulations clearly confirm that three phases coexist in rough pressure equilibrium, 
with temperature around 10 2 ,10 4 ,10 6 K, separately, and the hot gas contains a significant fraction 


of volume but little mass (de Avillez & Breitschwerdt 2004; Joung & Mac Low 2006 Joung et al. 




































-4- 


2009, Creasey et al. 2013; Gent et al. 2013 Walch et al. 2014). McKee & Ostriker (1977) applied 


their theory primarily to the ISM around us. We now want to extend to very different environ¬ 
ments, where the gas density and SN rate vary significantly from that in the solar neighbourhood, 
such as the Galactic halo, a starburst region, etc, where the multiphase structure of the ISM can 
vary vastly. In particular, when the SN rate is high enough, an equilibrium state no longer exists 


and wind generation is inevitable (e.g. Mac Low & McCray 1988; Scannapieco 2013 Gatto et al 


2015). 


An aspect of SN feedback that has long been neglected is that a significant fraction of SN 


progenitors are the so-called “OB runaways” (e.g. Gies & Bolton 1986 Stone 1991 Hoogerwerf 


et al. 2001 de Wit et al. 2005 Tetzlaff et al. 2011), thus SN can explode far from star formation 


sites. OB stars acquire this velocity when they reside in close binary systems and their companions 


explode as SN (Zwicky 1957, Blaauw 1961), and/or when they are in a crowded region of star 
clusters and get ejected due to dynamical interactions (P oveda et al.| 1967). Stone (1991) observed 
that 46% of O stars and 10% of B stars are in the “runaway” category with er ~ 30 krn/s. This 
means that OB stars can travel 50 —500 pc before exploding as SN. SN can therefore deposit energy 
outside dense molecular clouds, and may go into the low density inter-spiral arm or halo. The spatial 
distribution of the pulsars independently confirms this idea: the scale height of pulsars at birth is 


much larger than that of the Galactic HI layer, (cf. Narayan & Ostriker 1990). Consequently, the 


SN energy is much less prone to radiative losses due to low gas density. Type la SN, which are 
associated with the older stellar population with h ~ 300 pc, have similar effects. The “runaway” 
OB stars together with Type la SN are thus very likely to contribute to the launching of galactic 
outflow out of proportion to their intrinsic frequency (Ceverino fe Klypin||2009 Kimm fe Cen|[2014 


Walch et al. 2014), but this effect has been explored very little. The velocities of OB stars also 


determine the spatial distribution of the core collapse SN, which is important for the resultant 


ISM properties. For example, Gatto et al. (2015) found that random positioning of SN leaves most 


volume hot and thermal runaway, while SN exploding at density peaks result in little hot gas. 

In this paper we explore the interaction between SN and multi-phase ISM. The primary ques¬ 
tions we wish to address in this paper are: 


a. Given a multiphase ISM with a significant volume hot, how does a SNR evolve and impart 
feedback? In particular, how does the feedback differ from that in a uniform medium? 

b. Given a mean gas density n and SN rate density S, what is the resultant ISM, especially 
/v, hot? 

We propose three sets of experiments: First, as a baseline, we evolve a single SN in a uniform 
medium, to compare with the previous studies and parametrize the feedback effects. Second, 
we study how a SN evolves and imparts momentum and energy in an inhomogeneous medium, 
and compare it with the uniform case. Third, we allow multiple SN to shape the medium. We 
investigate the parameter space (h, S ) and analyze the resultant multiphase structure of the ISM. 
To be self-consistent, the multiphase ISM models in the second experiment are taken from the 
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third. 

The paper is organized as follows: In Section 2, we describe the code we use and relevant 
physics included; in Section 3, we present the single SN feedback in uniform and multiphase media; 
in Section 4, we show how multiple SN shape the multiphase ISM; we discuss our results in Section 
5, and conclude in Section 6 . 


2. Numerical Methods 


The numerical experiments are performed on a uniform-mesh code which uses the shock¬ 
capturing total variation diminishing (TVD) scheme (Ryu et al.| |l993) to solve the hydrodynamics 
equations. The ideal gas law is applied with the adiabatic index 7 = 5/3. Cooling and heating 
are calculated using the Grackle library^] as described in Smith et al. (2008/2011). Metallicity is 
assumed to be solar, that is, 2% of gas mass is in metals. Cooling and heating rates due to H, He 
and metals (including fine-structure lines) at a given density and temperature are pre-calculated 


using the Cloudy code (Ferland et al. 1998), which determines the ionization equilibrium solution 


for atom/ion species computed under an incident spectrum. In reality, the radiation held in a 
star formation region is very complicated and time-variable, and to be self-consistent one needs to 
appeal to radiative transfer. For simplicity, we do not take into account the ionizing radiation for 
this study. Photoelectric heating (PEH) through dust is included. We do not consider the heating 
of neutral gas due to cosmic rays or X-rays, since it is an order of magnitude smaller than the 


PEH (Draine 2011). At each time step (determined by the Courant condition), Grackle updates 


the internal energy of each cell by integrating the tabulated cooling/heating rates over that time 
step. Self-gravity, thermal conduction, molecule formation and cooling, and magnetic fields are not 
included in the current simulations. 

We test the code against the classic Sedov-Taylor solution with cooling turned off and a reso¬ 
lution of 1 pc. The overall agreement is satisfactory: the shock radius agrees to within 10% of the 
analytical result (Taylor 1950; Sedov 1959). 


3. Single SN evolution and feedback 
3.1. In uniform media 

In this section, we present the simulations of a SNR in a uniform medium, and analyze how a 
SN imparts its energy and momentum to the gas. We also compare our results with other recent 
numerical works. 


1 http: //grackle.readthedocs.org/en/latest/ 


















3.1.1. Fiducial run: Setup and Results 


We set up our fiducial run as follows: initially the gas is uniform with n = 1 cm~ 3 and 
T = 8 x 10 3 K. We inject E$n = 10 51 erg as purely thermal energy and 10 Mq mass into a sphere 
which encloses about 40 Mq of the ISM. The initial sphere has a radius of 7.5 pc. Both energy 
and mass are evenly distributed. The spatial resolution is 0.75 pc, and the box is 225 pc on one 
side. Periodic boundary conditions are applied, although the SNR does not reach the boundary. 
Background PEH is set to be 1.4 x 10 -26 erg/s per hydrogen atom, which is similar to that in the 
solar neighbourhood (Draine 2011). 


Fig. [ljshows the physical quantities (spherically-averaged) as a function of radius. Dashed lines 
in the upper-left panel indicate the shock radius as predicted by the Sedov-Taylor (ST) solution 
for the 7 timesteps we show. The position of the shocks in our simulation coincide with theory 
very well in the energy conserving phase. The peak density at this stage is ~ 2.5 n, somewhat 
smaller than that in the ST solution, i.e. 4n. This is likely due to the finite resolution and the 
application of spherical-averaging. After cooling becomes important, the blast wave slows down 
and lags behind the ST solution. A thin, dense and cool (T 10 4 K) shell forms, and the SNR 
enters the “pressure-driven snowplow” phase. We define i?o.8Eth as the radial position of the shock 
when 20% of the thermal energy in ST stage has been lost due to cooling. We find Ro.SEth = 22.1 pc 
at t = 4.1 x 10 4 year, in good agreement with the analytical calculation by 


Draine (2011 


(1) 


tc ool = 4.9 x 10 4 year( : 


n 


\—0.55/ 


Estss 


\ 0.22 


(2) 


' 1 cm 3 10 51 erg 

The equations are obtained by assuming the cooling rate A oc T~°- 7 for 10 5_7 3 K gas. Note that in 


Draine (2011), R C ooi and f coo i are defined somewhat differently, as the radius and time when 30% 


of the SN energy has been radiated away. 


While the shell has cooled and slowed down, the inner bubble is still filled with gas that is hot 
(T > 10 6 K) and fast-moving (v > 100 km/s). During this stage, the density of the hot bubble 
drops rapidly. This is because the fast-moving gas keeps running into the cool shell where it radiates 
energy very efficiently, which leads to substantial mass loss from the hot bubble. Accompanying 
the mass loss is a sharp pressure decrease, to even below that of the ambient gas. The precipitous 
decline of density and pressure makes the central engine quickly run out of steam. Thus, after a 
short-lived “pressure-driven snowplow” phase, a SNR enters the “momentum-conserving” stage. 
Without a propelling force, the cool shell moves forward due to its own inertia. The pressure of 
the shell, which is now beyond both the ISM and the inner bubble, drives itself to become thicker. 
Now the thick “shell” dominates the evolution of the SNR: while the outer edge of the shell still 
expands with a forward shock, the hot bubble it encompasses shrinks in size simultaneously. As 
shown in the temperature panel (lower-left), the radius of the hot bubble (T > 2 x 10 5 K) barely 
reaches 42 pc, while most swept-up mass moves to R > 80 pc. We define <fi as the ratio between 
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the maximum size of the hot bubble and R c 00 i, be. 


0 — -^hot ,max/Rcool • 


We find cj) = 1.77 for the fiducial run. The parameter 0 is important for the later calculation of 
the hot gas filling factor /y^ot (see Section [4]). Eventually, the velocity of the front shock drops 
to the sound speed of the ISM, and the temperature of the shell decreases until the PEH balances 
the cooling. The bubble created by SN will be filled with the ambient gas again, as seen from the 
increase of the density of the central bubble at the end of the simulation. 


Regarding the feedback effects, we are interested in how fast the blast wave propagates in the 
ISM, and how much energy and momentum it injects into the ISM as the shock reaches different 
radii. We present the results in Figj2j The radius of the blast wave R p , defined as the radius 
where the density peaks, follows a piecewise power-law function with time. The turnover happens 
at R p ~ Ro.sEth- We define the power-law indices for the ST and the post-ST stage as t]st and 
rj p_sT, respectively, i.e. 

R p oc t VST , for R p < Ro.sEth, (4) 

R p oc f ?p - ST , for R p ^ Ro.sEth- (5) 


A simple linear fit on log-log scale using the ordinary least squares method yields: t]st = 0.445 ± 
0.020 and ?? p -st = 0.292 ± 0.004. The ST solution predicts 77 st = 0.40. Chevalier (1974) found 
Vp-ST ~ 0.31; Blondin et al. (1998) had a somewhat higher value 0.33. 

The plot of energy vs R p (middle panel of Fig. [ 2 ]) indicates that the decline of kinetic energy 
Ek happens after that of thermal energy Eth• We use the parameter /3 to characterize this delay: 


@ — Ro. 8 Ek/Ro. 8 Eth, 


( 6 ) 


where Ro.sEk is defined similarly as Ro.sEth- Our simulation indicates Ro.sEth = 22.1 pc, Ro.sEk = 
30.5 pc and thus f3 = 1.38. The numerical value of (3 is understandable, since this process can be 
seen as an inelastic collision between the cool shell and the ISM, during which Ek only undergoes 
significant loss when the sweep-up mass is similar to that of the shell (38% increase in radius means 
roughly a 162% increase in volume and thus, mass). The decline of E to t, Eth and Ek seem to occur 
in a power-law fashion, and we define the corresponding indices a as follows: 


E tot « Rp Etot , 

for R p > Ro.sEth, 

(7) 

E t h oc Rp Eth , 

for R p > Ro.sEth, 

( 8 ) 

E k oc Rp Ek , 

for R p > R 0 . 8 Ek- 

(9) 


We find «Etot = —3.30 ±0.04, ctEth = —7.52 ±0.21, ctEk = —3.16 ±0.03. The loss of thermal energy 
is much more rapid than kinetic, indicating that at later stages of SN evolution, kinetic feedback 


dominates over thermal feedback (e.g. Thornton et al. 1998, Walch & Naab 2015 Simpson et al. 


2015. Kim Sz Ostriker 2015). We note that «Ek ~ —3 implies momentum V = (2mE ^.) 0,5 oc 






















(R 3 p E k ) 0 - 5 ~ constant, that is, momentum conservation starts at R p > Ro.sEk ~ 30 pc. This is 
confirmed by the V vs R v relation shown in the right panel of Fig. [2j 

The momentum V shows a power-law rising with R p , followed by a brief turn-over, and then 
reaches a plateau when R p > i?o.8Ek- Again, we define the power-law index yp for the rising phase 
of the momentum as: 

V (x R? p v , for R p > R 0 .8Ek- (10) 

We find yp = 1.21 ± 0.03. For the ST phase, when E k is constant and m oc R 3 , a simple scaling 
relation shows V = (2 mE k ) 0 ' 5 oc R p 5 . The final momentum Rend is 5.4 x 10 43 erg-cm/s, consistent 
with Kim & Ostriker (2015, hereafter K015). 


3.1.2. Dependence of results on input parameters 

The above scaling relations, especially the indices a, (3, y and numerical values of V e ndi are 
very important for SN feedback models. We want to test how the results change as we vary the 
numerical resolution, as well as with other model inputs, since in real galaxies, SN explode in 
various environments with different ISM density, PEH rate, metallicity, dumpiness, etc. In this 
subsection we study the effect of gas density, PEH rate and numerical resolution, and leave the 
task of exploring ISM dumpiness to later sections. 

In Table 1 we listed the models and the outputs. For n = 1cm -3 , we run three additional 
simulations with resolutions of 0.25 pc, 0.46 pc and 1.5 pc, respectively. The output parameters 
show variance within 10 %, and there seems no systematic change of the outputs as we improve the 
resolution by a factor of 6 , with the exception that ccEth seems somewhat higher in lower resolution 
runs. For the two highest-resolution runs, we have only lower limits on f> due to the limitation on 
the box size. 

Comparing four runs with different gas densities, i.e. n = 0.005, 0.1, 1.0 and 10.0 cm -3 , we 
find that a denser ISM shows somewhat more rapid decrease of E t ^ with R p , due to the more 
efficient cooling. But overall the difference is pretty small: decreasing n from 10 to 0.1 cm -3 makes 
cnEth drop by only about 20%. KOI5 also mentioned that the evolution is basically identical for 
different values of n when the radius and time are scaled relative to their values at shell formation. 
For the n = 0.1 and 0.005 cm -3 models, r] p sT (0.33) is larger than in higher density runs, which 
usually have 0.28 — 0.30. </>, which quantifies the volume of hot gas created by one SN, is in the 
range of 1.5-1. 8 . 

Vend is very insensitive to n. As presented in K015, this can be understood using a simple 
analytical argument. Since 


Vend = (2mE k )°- :> = [ 87 t /3 nm, H (l3Rcooi) 3 E k } 0 - 5 


( 11 ) 
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where E\ is simply a constant before momentum-conserving phase. Combining Eq. [T] and |TT| yields 

yt3,. ____ /_ / 11 0.13/ ^SN 


'Pi end ~5x 10 43 g • cm/s ( 


1 cm 


-3 


'-0.13/ xQ.94/ \1.5 

> 4o 51 erg ; ' IA 


( 12 ) 


In other words, in a uniform medium, SN (assuming a fixed total energy) inject an almost constant 
amount of momentum, independent of the density of the ISM. Our results agree with the analytical 
calculation fairly well. ctEk is somewhat larger than -3 for n = 0.1 cm -3 , which implies that 
momentum is not quite conserved at this stage, but is still slowly increasing. In this case the 
hot inner bubble has some pressure left to push the shell. V cn d is measured when the post-shock 
velocity drops to the sound speed of the ambient gas, 10 km/s. 

We test how the PEH rate affects the result for n = 10 cm -3 . From the data, we find there is 
no significant difference in the output parameters. One may notice that when we adopt a PEH rate 
1000 times the MW value, «Eth is somewhat smaller, which may seem unintuitive. This is because 
when we calculate £?th of the blast wave, we deduct the ISM energy (calculated by taking a sample 
outside the shocked region) from the total thermal energy enclosed in the blast wave domain. The 
high PEH rate keeps the ISM at ~ 10 4 K, while for the MW value, heating cannot balance cooling 
and the temperature of the ISM decreases with time. This makes the remaining E t ^ inside the 
blast wave smaller for higher PEH rates. 

Compared with K015, where they use a resolution of 0.25 pc for n = lcm -3 , we find very 
good agreement for «Ek 5 %>-ST and Vend- They have a somewhat smaller «Eth ~ —5.7, which can 
be due to the code and/or cooling curve we each adopt (they approximate the cooling curve by 
piece-wise power-law functions). 

In short, the evolution and feedback of a SNR in a uniform medium, as a function of time or 
R p , seem to be well characterized by (piece-wise) power-law functions. Overall, the dimensionless 
numbers a, /?, r/p-st and 0 vary little over a large dynamic range of n. 


3.2. In hot-dominated multiphase media 

In this section, we study how a hot-dominated multiphase medium (HDMM) affects SN feed¬ 
back. For a medium that is predominantly two-phase (warm/cold), one can refer to the recent 
works by K015. 

In principle, one needs many parameters to fully describe a multi-phase medium, for example, 
the size, density, shape of a single cloud; the location, number density and topology of clouds; 
the pressure of the inter-cloud medium; turbulence/velocity fields, etc, and to different extent, 
they all impact the evolution of a SNR. To exhaust the parameter space is not feasible. On the 
other hand, the above properties are not independent of each other, and more importantly, they 
are tightly correlated with SN feedback itself. (Other conditions, such as gravity, stellar outflows 
and radiation feedback, contribute to the ISM conditions as well, but arguably not as significantly 




-li¬ 


as SN, and to consider these effects is beyond the scope of this paper.) For self-consistency, the 
ISM models adopted in this section are based on the results from Section [4j where we explore how 
multiple SN collectively shapes the ISM. Other recent works use different methods to construct 
multiphase media. For example, K015 establish the warm/cold media through thermal instability; 


Martizzi et al. (2015) impose a log-normal density structure; Walch & Naab (2015) use fractal 


molecular clouds, with some runs having pre-ionization. 


3.2.1. Model Description 

We evolve a SNR in four inhomogeneous ISM set-ups, which have the same mean density 
1 cm -3 but different volume fractions of hot gas /y ,hot , including two idealized medium and two 
self-consistently generated ISM by multiple SN: 

i. “fvhO.OOa” (“fvhO.60” denotes /y^ot = 0.6, and “a” artificial) ; 

ii. “fvh0.90a”; 

iii. “fvh0.66ss” (“ss” denotes self-consistent); 

iv. “fvh0.82ss”. 

In the artificial ISM models “fvh0.60a” and “fvh0.90a”, clouds are identical to each other, and 
each cloud is spherical with a core of cold medium and a layer of warm phase. All three phases are 
in pressure equilibrium. The warm and cold phases are thermally stable. The size and mass of the 
clouds in “fvh0.60a” are chosen based on the results of the run nl_S'200 from Section [4] whereas 
“fvh0.90a” represents a more extreme model where 90% of the volume is hot and 90% of the mass 
cold. The locations of clouds are random. Detailed parameters of the three phases, such as the 
cloud radius R , c \, density n, pressure P/ks and volume/mass fraction fy/' f m are listed in Table [2] 
(a). The box size is 300 pc, and the resolution is 0.5 pc. The medium is static initially. 

For the self-consistently generated ISM, “fvh0.66ss” is the output from the run nl_S'200 of 
Section |4] at t = 1.1 x 10 8 year, and “fvh0.82ss” from nl_S'3000 at t = 7.3 x 10 6 year. Statistics of 
the ISM such as the volume-weighted density ny, pressure Py/ks and fy / f m of each phase are 
shown in Table [ 2 ] (b). To avoid the new SNR propagating outside of the simulation domain, we 
duplicate eight times the periodic box, assembling the eight identical boxes to make a larger one, 
which is 286 pc on each side, twice the size of the original one. The resolution 1.12 pc are the same 
as the original runs. 

We start by injecting 10 M 0 mass and 10 51 erg energy into a sphere with Mi n ; t ~ 200 M 0 
(Rinit ~ 12.5 pc) at the center of the box. For the artificial medium, the initial energy partition is 
30% kinetic plus 70% thermal. The injected thermal energy and mass are distributed evenly within 
the sphere. The initial E k of each grid is proportional to r 2 , where r is the distance of grid from the 
explosion center. We also tried cases where we inject energy as pure thermal or pure kinetic, and 
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the results converge after the blast wave reaches about 25 pc. For a self-consistently generated ISM, 
the energy added is purely thermal, and the injected mass and energy are uniformly distributed in 
the initial sphere. 


3.2.2. SN feedback in a medium with /y^ot = 0.6 

We first take a detailed look at the SNR in the “fvh0.60a” model, to understand the physics 
of the evolution in a multiphase medium, and then compare results of all models in this section. 

Fig. [3]shows 2D snapshots centered on the SNR in the “fvh0.60a” model at t = 1.9 x 10 4 ,1.1 x 
10 5 and 2.6 x 10 5 years, respectively. The first two rows use zoomed-in views. Visually the evolution 
is very different from that in the uniform ISM model. At t = 1.9 x 10 4 year, part of the blast wave 
has run into warm clouds where the shock velocity slows down. The geometry of the remnant 
thus deviates from spherical symmetry. The blast wave continues to travel fast in the hot inter¬ 
cloud medium - the least obstructed channel. At the same time, the shock heats and strips the 
clouds. The warm clouds are relatively easy to disrupt and accelerate, and are mixed with the hot 
component. The dense cold clouds remain largely intact and static due to their larger inertia and 
smaller cross sections. There is a visible density enhancement at the shock front in the hot gas, but 
a dense and cool (~ 10 4 K) shell, which features the “pressure-driven snowplow” phase of a SNR 
in a uniform medium, never forms. The post-shock region is also dissimilar to the Sedov-Taylor 
solution: the physical quantities inside the blast wave domain (excluding the clouds) are nearly 
uniform - it does not show the typical radial velocity distribution as v oc R, nor the progressively 
lower density and higher temperature toward the center of the blast wave. The shocked regions are 
turbulent, as a result of the violent interaction between the blast wave and the embedded clouds. 
At the end of the simulation, the warm clouds within R ~ 50 pc are largely disrupted, while the 
cold cores have survived. Clouds in the outer part of the SNR are less distorted, as the power of 
the blast wave subsides when the energy is spread into a larger volume. 

Fig. 0 shows the radial position of the shock R p vs time and the energy/momentum feedback 
as a function of time and R p . For comparison, the dashed lines indicate the SNR evolution in a 
uniform model with the same mean density n = 1cm -3 , same as what is shown in Fig. i It is 
obvious that the evolution in a HDMM is very different from that in a uniform case. 

Panel (a) of Fig. [4] shows the evolution of R p vs t. We determine R p in the following way: 
since initially the medium is static, we divide the simulation domain into concentric spherical shells 
centered on the SN explosion site with the thickness of the cell size, and define R p as the radius 
beyond which no shells has more than 10% of cells with v > 10 km/s. From the plot, the blast wave 
travels much faster and further in the multiphase medium. For the uniform case, the power-law 
index of R v -t turns from 0.4 to ~ 0.29 after the SNR enters the “pressure-driven snowplow” phase; 
in contrast, in the multiphase model is more close to a single power-law function, with a slope ~ 
0.55. As a result, by t = 1.8 x 10° year, the difference in shock radius is more than a factor of 3, 
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and the volume enclosed differs by a factor of ~ 30. The hot-dominated medium transfers the blast 
wave out much more efficiently. 

To see the importance of the clouds in determining the blast wave velocity, we plotted R p —t for 
a uniform medium with n = rthot = 4 x 10“ 3 cm -3 and n = n warn i = 0.4 cm -3 , respectively (dotted 
lines in panel(a) of Fig. [4]), as predicted by the ST solution up to their separate t coo i (Eq. [2j. At 
t < 10 5 year, the curve for “fvhO.OOa” lies in between the above cases; after that, the evolution 
almost coincides with that in a n = 4 x 10 -3 cm -3 medium without energy loss. This is because 
that during the early stages, the blast wave is powerful and interacts more with clouds, thus the 
velocity is reduced, whereas at later stages, the shock largely bypass the clouds and travels more 
freely in the hot medium. This agrees with what we have seen in the slice snapshots. 


We estimate what fraction of a cloud interacts with the blast wave. For a cloud of the size 
2R C \, density p c \ and pressure P c \ that is located at a distance d c \ from the center of a spherical blast 
wave, at a time t, the length scale of the cloud passed by the shock is R s h(t) = f to ( R =d ^ v s h(t')dt ', 
where the shock velocity in the cloud (t) V P sh(,t)/Pd ~ Cs tC iy/P sh (t) / P ch P sh is downstream 
pressure of the shock, and c SiC i the sound speed of the pre-shock cloud (Klein et al. 1994). Assuming 
P s h is uniform within the SNR and oc P~ 3 , and R p oc t 71 , then at t = tf a d e , when the blast wave 
decays to a sound wave, i.e. when P s h ~ P c i, we have 


Psh 


2 P c i 2 — 3r/ t s , cl 


( ^ fade ) [1 _ ( ^ cl ^ (2—3r))/2r)j ^ 


Pfade 


(13) 


where P s h/2P c i is the fraction of a cloud being shocked, f SiC i is the sound crossing time of the 
cloud (= 2P c i/c Si d), and Pf a de is the radius of the blast wave at tf a de- Note that the closer a cloud 
is to the center of the SNR, the larger R s h/2R C \ it has. Given the approximate value i] ~ 0.55, 
Pfade ~ 200 pc, we have P s h/2P c i ~ 0.5 and 1.2 for the cold and warm cloud, respectively, at 
d c \ = 0. This means that only the warm clouds near the center can fully interact with the blast 
wave, and cold clouds are mostly bypassed. In other words, for most cold clouds encompassed by 
a SNR, the interior of the clump does not even know the existence of the blast wave before it fades 
away. 

Panel (b) of Fig. [4] shows the SN energy as a function of time. The SN energy is obtained 
by summing up all energy within R p with the energy from the enclosed ISM deducted. The ISM 
energy density is acquired by averaging a sample outside of the SNR. The evolution of energy has 
two stages: an energy conserving phase followed by a radiative phase, divided by when roughly 20% 
of total energy is lost. For “fvh0.60a”, this transition happens at t ~ 5 x 10 4 year and R p ~ 65 pc. 
The time is similar to that in the uniform case, while the corresponding R p is about 2.7 times larger. 
The energy feedback shows interesting contrasts to the uniform case for both stages. Initially, after 
the deposition of energy as 30% kinetic and 70% thermal, Ek rises while Pth declines, deviating 
from the ST solution. E k reaches its maximum of about 40% of SN energy at approximately the 
end of the energy conserving phase. The results do not show qualitative variation when we change 
the initial energy deposition to 100% thermal or 100% kinetic - E k always exceeds that in the 
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ST solution, though Ek/E tot when 1% peaks can vary by about 20%. This initial conversion of 
thermal to kinetic energy is reminiscent of free expansion, and is most likely due to expansion of 
the shock-heated clouds into the tenuous medium. 

Despite that the decline of £th happens earlier than the uniform case, the rate at which it 
decreases is much smaller in a HDMM. While in the uniform medium, £% drops by a factor of 50 
within 2 X 10 5 years, the remnant in “fvhO.OOa” loses only about 50%. The main reason for the 
significantly reduced cooling rate is that the blast tends to choose a more tenuous path, which is 
cooling-inefficient, and largely bypassing the dense clouds. Another reason for the mild loss rate of 
E t h is that at later stages E^ can convert to £% through turbulence. E^, on the other hand, does 
not show very different time evolution from that in the uniform phase. 

The plots of energy vs R p (panel (b)-(d)) contain no new information but simply re-present 
the data from R p — t and E — t. It is obvious that the energy is transported to a much larger volume 
by the fast-moving blast wave, and drops much more mildly with R p in a HDMM, especially for 
E t h- In a uniform medium, SN energy is negligible after the catastrophic cooling at R p = 65 pc, 
whereas in a multiphase ISM, when the shock fades to a sound wave at R p > 130 pc, ~ 25% of the 
SN energy is still maintained in the “fvh0.60a” model. 

Momentum V as a function of R p is shown in panel (f) of Fig. [4j The momentum does 
not increase as fast with R p in a multiphase ISM. Given V = \J2mE^, the smaller V and larger 
E^ implies much less mass is involved in sharing the kinetic energy in a HDMM. Note that no 
thin shell forms at the shock front which carries most of the radial momentum, as happens in 
a uniform medium. At the end of the simulation, the radial momentum in “fvh0.60a” is about 
4.3 x 10 43 g-cm/s, and is still rising. But it is not likely to keep increasing for long, since the post¬ 
shock pressure is about to drop to that of the ISM and will no longer convert into the outgoing 
momentum. This suggests that the total momentum injected in a clumpy medium is close to that in 
the uniform case 5.3 x 10 43 g-cm/s. We also plot the total momentum (dotted line), which initially 
coincides with the radial momentum, but later exceeds it since a non-radial component develops 
as the blast wave interacts with the clouds. At the end of the simulation, the total momentum is 
roughly 30% larger than the radial one in “fvh0.60a”. 


3.2.3. SN feedback as a function of fy^ ot 

In this subsection, we compare SNR evolution in all 4 multiphase models, and present the 
fitting formulae to quantify SN feedback in a HDMM. The fitting formulae in this section are for 
F?SN = 10 51 erg, and gas average density n = 1 cm -3 and fy^ot > 0.5. 

Fig. [5] shows the snapshots of the three multiphase ISM models “fvh0.66ss”, “fvh0.90a” and 
“fvh0.82ss” at t = 3 x 10 4 year. It is clear that the blast wave travels mostly in the hot phase. For 
“fvh0.66ss” model, where the clouds form walls and hot bubbles do not fully connect to each other, 
the shock goes further in the low-density bubble and deviates from spherical symmetry. Within the 
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SNR, the pressure is higher toward the clouds. For “fvh0.90a” and “fvh0.82ss”, where the clouds 
are smaller with nearly random locations, the blast waves are close to spherical. 

Fig- © compares SN feedback in all multiphase ISM models. For reference, the uniform case is 
shown with the blue dashed lines. For “fvh0.90a”, R p , Eth , Ek are determined the same way as for 
“fvh0.60a”. For the self-consistently generated ISM models, which are turbulent and the SNR may 
have an irregular shape, we briefly describe how we determine the size and energy of the SNR. For 
“fvh0.66ss”, we adopt the following criterion to determine if a cell is within the SNR. A post-shock 
cell should satisfy at least one of the two conditions: a. v > 150 km/s; b. nT > 10 4 ' 5 cm~ 3 K. 
The adopted velocity and pressure floors are beyond (but very close to) the maximum values in the 
pre-SN ISM. We then define the “effective radius” R p = (SFshock/^vr) 1 / 3 , where F s hock is the total 
volume of the post-shock region. For the “fvh0.82ss” model, since the ISM has higher pressure and 
velocities, and the spatial variations are larger than “fvh0.66ss”, it is difficult to obtain R p in the 
same way. Taking advantage of the relatively regular shape of the SNR, we take snapshots and 
identify the R p by eye. The SN energy is determined by running a reference simulation for each 
self-consistent ISM model where no SN is added. At each time step, the difference of Ek and E t h 
between the two simulations are deemed those associated with the SN. 


In Fig. [6j one may first notice that the four multiphase models show roughly similar SNR 
evolution, but they are quite different from the uniform model (except panel (e) Ek-1 ). The compa¬ 
rable evolution is interesting because the four multiphase models have fairly different configurations: 
cloud shape, /v.hoti turbulence, and even resolution. This implies that those differences are sec¬ 
ondary to the fact that they share the same feature of a HDMM, distinct from the uniform case. 

Let us then look at the feedback in some details. Panel (a) of Fig. [6] shows R p vs t. SNRs 
in all multiphase ISM models go much faster and further than the uniform medium. For both 
the artificial and self-consistent ISM, R p is larger for higher /y^ot a t a given time, but overall the 
difference is small. For “fvh0.66ss”, at later stages, R p reaches its maximum of 95 pc, when the 
shock subsides and merges with the cloud boundaries of the bubbles. The SN feedback is confined 
by the walls, thus it does not influence as large volume as when hot gas has a connecting topology. 
Overall, R p vs t is roughly represented by the following power-law function 


R p = 22 pc( 


t 


\0.55 


7 x 10 3 year ; 


(14) 


The power-law index, 0.55, is larger than the ST solution. The corresponding curve is shown 
with the black dotted line (same linestyles are adopted for other fitting formulae presented in this 
subsection in Fig. [6]). 

Panel (c) of Fig. [6] shows the time evolution of E t h, which is characterized by an early onset of 
decline and a small cooling rate for all multiphase models. Self-consistent ISM models lose -Fth at 
a somewhat larger rate, which may be due to their turbulent nature that enhances the dissipation 
rate. If we compare the two artificial ISM models or the self-consistent ones, a larger /y^ot yields 
a smaller cooling rate. -Eth(t) has a power-law-like form, with power-law indices (absolute value) 
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much smaller than the uniform case. We adopt a very simple formula to fit the power-law indices 
for the two self-consistent models, which are inversely proportional to /v,hot : 


E t h (t) 


0.65£sn for t < 0.09t C ooi,fi) 

0 - 65S SN(wwJ-)—0.36x(0.8// v ,hot) for t > 0 .09t coo i ifi 

U.U9t coo in 


(15) 


where t coo i,n = 4.9 x 10 4 year(h/l cm~ 3 )~°- 55 (EsN/10 51 erg) 0 - 22 , following Eq. [2j Note that the 
functional form of the power-law index is chosen somewhat arbitrarily, with no direct physcial 
motivation. Only two values of /v,hot are fit, meaning that it is not strongly constrained; however, 
this is probably only of secondary importance for sub-grid modeling, as /y,hot does not vary strongly 
across HDMM conditions. We set a ceiling for it at 65% of the SN energy, to avoid the divergence 
at t = 0. The onset of cooling, which is about 10% of f CO oi,fn is much earlier than the uniform case, 
while the canonical value of the power-law index of E t ^ — t for a HDMM, -0.36, is significantly 
smaller than the uniform case, ~ — 2. 

Panel (d) of Fig. [6] shows Eth vs R P • The fitting formula can be obtained from combining Eq. 
[I4land[l5l 


Eth(Rp) = 


0.65Esn 

0.65E S n( 


Rit 


0.73 E, 


cool,n 


for R p < 0.73i? coo i )fi , 
0.65x(0.8// v ,hot) for Rp > 0.73E coo i, h , 


(16) 


where 72 coo i ifl = 23.7pc(h/l cm~ 3 ) _a42 (EsN/10 51 erg) 0 - 29 , following Eq. [I] Again, the SN energy 
starts to cool at a somewhat smaller R p , while the dependence of E t h on R p is much less sensitive 
in a HDMM than in a uniform medium, which has a power-law index of about -7. Note that the 
evolution of E t h in our three-phase ISM simulations is different from that in two-phase models. 

find that the the power-law slope of the declining phase of E t h vs R p is 
similar to the uniform case (see their Fig 4). K015 show that the energy evolution in a two-phase 
ISM is relatively close to the uniform case with the warm phase density (see their Fig 10), since the 
cold chimps only have a small volume fraction. In a three-phase medium, we find that evolution 
is more complicated. Although the late evolution of R p is close to that in the hot medium only, 
the energy evolution cannot be rescaled to the uniform case with any density. It thus suggests that 
some new model is needed which takes into account both hot gas and warm clouds. 


Martizzi et al. (2015 


Panel (e) of Fig. [6] shows the time evolution of E^ , which in multiphase media does not show 
large deviations from that in the uniform medium, although in the energy conserving phase, Ek 
has a higher fraction of SN energy than 30%. A simple fitting formula is 


Ek(t) 


0.35E S n 

0.35E sn (———) -0 - 87 

^cool,n 


fol t ^ ^cool,n? 
for t > ^cool,n* 


(17) 


Efc vs R p is shown in panel (f) of Fig [6j Similar to E t h , E^ is conveyed to a much larger 
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volume in a HDMM. For the fitting formula, combining Eq. 14 and 17, we get 


Ek(R P ) 


0.35Esn for R p < 2.7i? coo i,n, 

°' Msn( 9TF- )_1 ' 58 fOT R P > 2 - 7i? cool,n- 

Z. (I X C ool.n 


(18) 


Note that in the declining phase, the power-law index of R p , -1.58, is less steep than the uniform 
case, -3, consistent with that the momentum is not conserved, but still rising when Ek declines. 


The momentum evolution is shown in panel (b) of Fig [6} Both the total momentum and the 
radial component are presented. We only plot the artificial ISM cases, since in a turbulent medium 
it is hard to isolate the momentum directly related to the new SNR. The SNR in the “fvh0.90a” 
model has a consistently smaller momentum than that in “fvh0.60a”, at a given an R p , since less 
mass is involved in the blast wave in a more hot-dominated medium. As the shock goes to further 
distances, momentum is developed in non-radial directions and the post-shock region becomes more 
turbulent. For the most part of R p , the momentum in a multiphase medium is much smaller than 
the uniform case, but the terminal momenta are likely to be similar. For the radial momentum, we 
adopt the following fitting formula: a power-law increase with R p followed by a plateau as in the 
uniform case: 


R rad = 


= f 5 x 10 43 g • cm s -i (t:5 |e_)0-57(^)-o. 5 for Rp < 7.5R cool , fi (%i) 0 - 88 , 

for R p > 7.5R coolifi (%^) 0 ' 88 . 


5 x 10 43 g • cm s 1 


(19) 


The power-law indices of EgN and n follow those of the uniform case (Eq. 11). Note that while a 
single power-law describes the evolution in “fvh0.60a” almost perfectly, it does not agree well with 
“fvh0.90a”, which may be better fitted by a broken power-law. Here we simply make /y,hot a factor 
that only influences the zero point of curve, to reflect that less momentum is developed at a given 
R p in a medium with a larger fy^ot ■ 


Note that we have shown the uniform case with the same mean density of the HDMM models. 
Since the cold phase that contains most mass actually fills little of the volume, from a physical 
perspective, one could also compare to uniform models with the same median density. We do not 
show them here to avoid overcrowding the figures. However, as we have shown in Section 3.1 the 
SN feedback in a uniform medium is well described by power-law functions, and the power-law 
indices vary little for different gas densities. So as a rule of thumb, the feedback curves simply shift 
horizontally for different densities. 


In summary, the evolution of a SNR in a HDMM is characterized by a much larger impact 
volume, a reduced cooling rate, and a slowly developed momentum. When the medium is uniform 
in density, there is no preferential path and the shock encounters mass at the same rate in all 
directions; when the medium is clumpy, however, the shock tends to choose the more tenuous 
regions where it travels much faster and further. The connected rarefied channels give vent to 
the blast wave energy, which largely goes around the clouds at later stages. Consequently, not all 
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mass within R p interacts significantly with the blast wave. A HDMM with connecting tenuous 
tunnels therefore greatly facilitates SN feedback. In cosmological simulations which cannot resolve 
the multiphase, turbulent medium, it has been commonly assumed that SN explode in a uniform 
medium with the mean cell density. The cooling rate can thus be severely overestimated. A more 
realistic and self-consistent model should consider a SN-modified multiphase ISM. 


We can define the domain of influence of a SN in 4D spacetime: 


V 4 D = -^--Ro.2*o.2- 


( 20 ) 


If we define -R 0.2 and to .2 as the R p and time when 20% of SN is energy is left, then for a 10' 51 erg SN, 
R 0.2 = 33.8 pc and to .2 = 1-7 x 10 5 year in a uniform medium with 1cm -3 , and -R 0.2 = 190.0 pc and 
to .2 = 3.5 x 10 5 year in a multiphase medium with same mean density and /y hot = 80%, according 


to Eq. 15 - 18 The impact volume is therefore 178 times larger and time 2.1 times longer, making 


a 374 times bigger V 4 D hr a HDMM than in a uniform medium. 


In real galaxies, it is usually the collective effect of many SN that shape the ISM and launch 
galactic winds. A larger V 4 D for a SNR means that a new SN is more likely to explode within 
the Rid of an old one. Thus the chance for SNRs to overlap significantly increases, and energy 
from multiple SN are “grouped”, which may result in a big runaway explosion. It is therefore very 
interesting to investigate the question “Under what circumstances will a HDMM form?”, which we 
will do in the next section. 


4. Multiple SN shaping multiphase ISM 

In this section, we allow multiple SN to shape the ISM. We experiment with two parameters, 
the mean density n and the SN rate S, and analyze the resultant ISM. In particular, we want to 
find out when there will be a HDMM, which has the potential to drive a wind. 


4.1. Model Description 


The gas is initially uniform with a density n, temperature 3 x 10 3 K and zero velocity. Periodic 
boundary conditions are applied. SN are assumed to have random locations (we test how good this 


assumption is later in Section 5.1 taking into account runaway OB stars). We explore the parameter 
space of (h, S ) to study the resultant ISM systematically. Each simulation box contains 128 3 cells. 
The box size of each simulation varies with n : the length is chosen to be 6R coo i (Eq. [TJ. The 
adaptive box size allows us to achieve higher resolution for a denser ISM, thus better capturing the 
evolution of SNRs and avoiding overcooling problems, which are frequently encountered in coarse 
resolution simulations. 


Each supernova is injected as 10 51 erg thermal energy and 10 Mq mass, which are distributed 
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evenly within a sphere. We choose the initial radius such that the enclosed mass is close to but 
less than M- m it, and we adopt = 200 Mq for n < 1 cm' 3 and 100 Mq for n > lcm -3 . As for 
the duration of the simulation, we note that there is an intrinsic timescale after which ISM has 
been fully bombarded by random-located SNR, so that the initial uniform distribution of gas is 
destroyed entirely. This is roughly 


tb = 


S (2i? coo i) 3 


= 9.2 Myr ( 


-plf 11 \1.26 

1000 kpc -3 Myr -1 Mem -3 


( 21 ) 


We therefore choose the simulation time to be i s i m = 2.5 4, so at the end of each simulation ~ 
66 SN have exploded in the box. The time span corresponds to about 10% of the gas depletion 
time (= nmn/(S • 100Mq)) for n = lcm -3 runs, and 40% for n = 30cm -3 . We find that the ISM 
reaches steady state (if there is one) long before the end of the simulation. 


We assume PEH rate scale linearly with the S, and adopt the following formula for a uniform 
PEH rate per H atom (Draine 2011): 


r pe = 1.4 x 10 26 erg • s 1 


100 kpc 3 Myr 


r 1 ' 


( 22 ) 


We implicitly assume a constant dust/gas ratio, in line with a constant metallicity. Given a nominal 
dust photoabsorption cross section of FUV photons per H atom, 10 -21 cm 2 (Draine 2011), we find 
that most of the simulation boxes are in the optically thin regime, except for the n = 30 cm -3 runs 
where the total optical depth across the box is about 3.3. Therefore a linear scaling relation of 
PEH rate with S is, in general, not a bad approximation. 


4.2. Choice of Parameters 


The (h, S ) combinations encountered in different galaxy environments are not well-constrained, 
and are likely to be time-variable. Although on scales of kpc and above, we have the Kennicutt 
relation on the star formation (SF) rate and mean gas surface density (Kennicutt 1998) (which has 
intrinsic scatters itself), on a smaller scale, like that of our simulation boxes, recent observations 
have shown large deviations from it. For example, by doing stellar counts for the star forming region 


in the MW, Heiderman et al. (2010) found that the SF rate can be above the Kennicutt relation 
by at least an order of magnitude for S gas ~ lOOM 0 pc -2 (although this is for small individual 
star-forming regions, not for large scale ISM). Besides, the Kennicutt relation uses the observables 
- the column densities of gas and SF rate, instead of the volume densities, and the scale height 
of the gas for external galaxies are usually not well-known. Furthermore, Type la SN, which are 
associated with the old stellar populations, do not correlate with gas properties. Lastly, Type II 
SN, because of the non-negligible velocities of their progenitor OB stars, can migrate out of the SF 
sites. All these contribute to the scatter of n — S relation in real galaxies on different length scales. 

Given our incomplete knowledge and the potentially big scatter in the relation of h and S, we 
treat the two as free parameters, the relation between them being determined by physical modeling 
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that depends on self-gravity, which has been ignored in this work (see, for example, Kim et al 
(2011, 2013) for self-consistent determination of Xsfr as a function of S gas and E star for disk 


galaxies; for elliptical galaxies, SN are mostly Type la from old stellar populations, therefore such 
self-regulations do not exist.). In terms of choosing what parameter space to explore, we try to 
combine both observationally- and theoretically-motivated regimes. From the observational side, 
we cover the following regions and their neighbourhood - (1) MW disk-average (n ~ 1cm -3 and 
S « 200 kpc -3 Myr -1 ); (2) MW halo: based on vertical distribution of SN and halo gas; (3) Scaling 
the MW disk-average to “starburst” regime using the Schmidt relation S oc p/td yn oc n L5 . Beyond 
the above domains, we also expand the parameter space after we analyze the results of the existing 
runs, to include those that are theoretically interesting. For example, given the density of the ISM, 
there is a critical SN rate S/nt above which no steady state exists. We list in Table [3] the basic 
parameters for all simulations in this section. The symbolic name of a simulation, for example, 
h0.1_510, means n = 0.1cm -3 and S = lOkpc -3 Myr -1 . “Res” is short for “Resolution”. 

For the rest of this subsection, we explain how we estimate the n - S correlation in the vertical 
direction of the MW disk. Observation of the 21 cm line indicates that the neutral medium lies 
within a flat layer approximated by: 

«hi = 0.57cm -3 {0.7 exp[— (z/127pc) 2 ] + 0.19 exp[— (z/318pc) 2 ] + 0.11 exp(—|z|/413pc)}, (23) 


where z is the vertical distance from mid-plane (Dickey & Lockman 1990). Gaensler et al. (2008) 


deduced from pulsar dispersion measure and diffuse Ha emission that the ionized warm medium of 
the MW follows an exponential function: 


n e 


= 0.014cm 3 exp(—|z|/1830pc). 


(24) 


The hot ionized gas is less constrained, and the density is very low, so we leave it out for this 
study. The sum of the above phases then yield the vertical distribution of the MW gas. 


The SN frequency for the MW is about 1/60 year 1 for core collapse SN and 1/250 year 


„-i 


for Type la (Cappellaro et al. 1997). Type la SN come from the old stellar disk which has a scale 
height of ~ 325 pc, therefore we adopt the following vertical distribution of Type la SN, as averaged 
over the MW disk (assuming 10 kpc radii in size): 


Si a = 19.6kpc 3 Myr 1 exp(—|z|/325pc). 


-l. 


(25) 


For core collapse SN, the vertical distribution is inferred from that of pulsars at birth, which is a 
sum of two Gaussian functions (Narayan & Ostriker 1990): 


See = 167kpc 3 Myr 1 { 0.75 exp(z/120pc) 2 + 0.25 exp(z/360pc) 2 }. 


-l 


(26) 


Putting Eq. |23||26| together, we get an n — S relation along the vertical direction of the MW. This 
is shown by the blue solid line in Fig. m and |13| that we will describe later. 



















Table 3: Input/Output parameters of simulations in Section 



TR = Thermal Runaway, Y = Yes, N = No, m = metastable state 
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4.3. Results 


4-3.1. Slices of some output ISM 

Fig. [7] shows slices of density, temperature, pressure and velocity of three simulation examples. 
The upper panels are for the run nl_S'200, corresponding approximately to the solar neighbourhood. 
Gas settles into three distinct phases with contrasting density and temperature: Hot T 10 6 - 6 K, 
warm T ~ 10 4 K, and cold T ~ a few hundred K, which have been observed by a variety of 
techniques and summarized in McKee & Ostriker (1977). These phases are in rough pressure 


equilibrium, with variations of a factor of ~ 30, significantly smaller than the dynamical range 
of density and temperature. Velocity fields indicate that the ISM is turbulent, with the lighter 
component moving faster. Hot gas fills about 60% of the volume after reaching steady state. These 


results agree with the previous analytical and and numerical simulations (e.g. McKee Sz Ostriker 


1977 

Joung & Mac Low 

2006; 

Gent et al. 

2013 

Gatto et al. 

2015) 


Middle panels are for n0.3_ST00, an example of a HDMM. SNRs overlap with each each other, 
rendering a connecting topology: ~ 90% of the volume is filled with hot gas, and dense clouds are 
bathed in it. The clouds mainly consist of the warm phase. The deficiency of the cold medium is 
attributed to the less efficient cooling of less dense ISM, and the PEH is able to keep the temperature 
of most clouds to rs_/ 10 4 K. 

Lower panels show the results of hl0_S'3000, corresponding to a star forming region. The SN 
rate is insufficient to retain a hot phase. Hot bubbles close up before they overlap with each other. 
SN do not play a significant role in shaping the thermal state of the ISM. The two-phase medium 
is maintained by the balance between cooling and PEH. On the other hand, the momentum input 
from SN drives turbulent motions in the gas, and is therefore important for shaping the dynamical 
state of the ISM, which may have strong implications for star formation (Mac Low & Klessen 2004 


McKee & Ostriker 2007). 


4-3.2. Evolution of the different phases 

(1) hl_5200 (MW disk-average): example of reaching steady state and pressure equilibrium 

Fig. H presents the time evolution of the three phases for the run nl_S'200. We show in the 
four panels the volume-weighted pressure and Mach number, and the volume and mass fraction 
of each phase. Each data point represents one snapshot. The properties of the hot gas exhibit 
larger fluctuation, reflecting periodic explosions of SN. The evolution of the warm and cold phases, 
on the other hand, are smoother. Initially, the gas is uniformly 3000 K; within a few Myrs, cold 
clouds quickly build up and occupy > 80% of the mass, and the cold and warm media both settle 
into their separate thermally-stable temperatures. Thereafter, the mass fraction and temperature 
for these two phases keep almost constant. The separation of the two phases happens as a result 
of the thermal instability, facilitated by perturbations from SN. The Mach number is higher for 
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colder phases: the cold and warm media are mildly supersonic with Mach number ~ 2.3 and 1.3, 
respectively, while the hot gas is subsonic (the slightly supersonic cold/warm phase and subsonic 
hot phase turn out to be the case for most of the runs, see Table [3]). All three phases are in 
approximate pressure equilibrium, with nT ~ 3 — 4 x 10 3 cm ~ 3 K. These results agree very well 
with the observation (e.g. Jenkins & Tripp 2011). Hot, warm and cold gas occupy 60%, 40%, 1.5% 


of the volume, respectively, consistent with the observations in the solar neighbourhood (e.g. Heiles 


& Troland 2003; Murray et al. 2015). Note that observations of 21 cm absorption-line yield about 


60% of the mass of neutral hydrogen is warm (Heiles & Troland 2003), which seems inconsistent 


with our result. However, we use T < 10 3 K as the criterion for being “cold”, different from the 
the observation, which is roughly T < 200 K. Also we do not distinguish molecular gas from the 
cold medium (nor do we track the ionization state) in our simulations. The mass ratio between 
cold HI and molecular is about 1.38 for the MW (DrainepOll). Taking this into consideration, the 
revised mass ratio of cold to warm in our simulations becomes 3.3 (even assuming all warm gas is 
neutral), still larger than the observed 0.67. Some effects that we do not include may account for 
the over-condensation of the cold phase, for example, UV from massive stars can ionize and heat 
the cold medium; thermal conduction by hot gas may evaporate some of the cold gas and form a 
warm component. Another possibility is that if a few SN happen to have exploded in the dense 
clumps, they can effectively turn the cold medium into warm and hot (e.g. Gatto et al. 2015). The 
stellar winds from their progenitors can have a similar effect. We will discuss more on the position 
of SN in Section 5.1. 

(2) n0.1_S'50 (MW disk-halo interface): example of a HDMM and “thermal runaway” 

The (n, S ) of this run corresponds to the disk-halo interface of the MW (see Fig. [IT] below). In 
the left panel of Fig 9, we show the volume-weighted pressure evolution for the three phases. Same 
as in Fig. 8 , each data point is calculated for one snapshot. Hot gas is over-pressured compared 
to the other two phases, which reflects what we have found in the previous experiment: the blast 
wave travels preferentially in the hot medium. It would take roughly the sound-crossing time for 
the clouds to re-establish the pressure equilibrium. The overlapping SN determine the overall 
properties of the ISM. The medium does not reach a steady state at the end of the simulation, and 
it is not obvious that it will ever reach one. This is because, despite the rising density of the clouds 
which enhances cooling, the volume is increasingly occupied by low-density hot gas where cooling 
is negligible; since blast waves lose the majority of their energy in the clouds, cooling is inevitably 
delayed until the shock sweeps enough volume. Similar phenomena have also been observed in 


some of the simulations in Gatto et al. (2015). Scannapieco et al. ( 2012 ) discovered a transitional 
ID velocity dispersion ~ 35 km/s for the turbulent ISM, above which the rarefied gas cannot cool 
efficiently and the ISM undergoes a thermal runaway. We can understand the thermal runaway 
state in this way: if the SN rate is high enough for multiple SN to overlap, the ISM will be multiphse 
where the hot bubbles connect to each other. As we have found in the experiment in which a single 
SN propagates in a HDMM, the volume encompassed by the blast wave is much larger, and the 
cooling is suppressed, thus the spatial-temporal domain of a single SN (V 413 ) is greatly increased. 
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This makes the overlapping of the SN reinforces itself, which leads the ISM to a “runaway” state. 
In real galaxies, the persistent increase in pressure implies that the gas will expand, and hot gas 
can blow out and form winds (which are not captured in the periodic box simulations). 

(3) hl0_S'3000 (SF region): example of a subordinate hot phase 

As shown in the right panel of Fig. [9j the ISM in this case has only warm and cold phases 
most of the time. Hot gas is present only briefly after SN explosions. SN are not determining the 
general thermal state of the ISM, since the pressure of the hot gas is significantly smaller than 
the other phases, and the time-averaged volume fraction is no more than a few percent (see Fig. 


“pressure-driven snowplow” phase, hot gas moves much faster than the shell, and quickly runs 
into it and loses energy, which leaves the hot bubbles at an extremely low density, and therefore 
low pressure (Fig. [I]). Note that right after SN explosions, the hot gas is over-pressured, but it 
only lasts till the shock becomes radiative, and the time interval is very short compared to our 
sampling. In fact, according to Eq. [ 2 J t coo i = 0.014 Myr for n = 10 cm” 3 , and yet the time interval 
between two consecutive SN is 2.08 Myr, and our sampling will not catch the over-pressured stage. 
In other words, if we could sample the data at a frequency larger than every 10 4 year, we would 
see the very spiky evolution of the hot gas pressure; if we could enlarge our box size (but keep 
S unchanged) so that at any given time, there are many SN exemplifying different stages of their 
evolution, we would see the mean pressure of hot gas go beyond that of warm/cold, because the 
spatially-averaging is biased towards regions with pressure that is orders-of-magnitude higher than 
the median. In other words, when we say hot gas pressure is lower than that of the other phases, 
we mean that at any given time, if we randomly choose a location that is above 2 x 10 5 K, it is 
most likely to be under-pressured. 


11 below). The reason for the under-pressure of the hot phase is this: once the SNR enters the 


(4) /v, hot ~ 0.6 ±0.1 : metastable state of the ISM and transition to “thermal runaway” 

For the simulations which reach a steady state and have /y,hot ~ 0.5 — 0.7, we sometimes 
observe an interesting transition from “steady state” to “thermal runaway”. Fig. [T0| exemplifies 
such a transition. This is for (n,5) = (1,200), with the same set-ups as the fiducial run. The 
locations of the SN are random, but different from the fiducial one. Initially the ISM reaches a 
steady state (same as Fig. [8]): Each of the three phases occupies a stable portion of the volume, 
and they are roughly in pressure equilibrium. Starting from t ~ 55 Myr, however, hot gas quickly 
comes to dominate the space, while the volume fraction of the warm phase drops from ~ 50% to a 
few percent. Meanwhile, the gas pressure rises progressively, and the equilibrium is disrupted - the 
hotter phase has a higher pressure. These phenomena are the same as those in a “thermal runaway” 
state. The transitions are somewhat stochastic, and can happen at any time after the steady state 
is established (the example shown in Fig. [8] is stable to the end of the run t = 110 Myr). The 
stochastic transition is not so surprising: the configuration of the ISM with half the volume hot 
is on the verge of changing from “isolated bubbles” to a “connecting” topology; any new SN may 
break the balance and lead to the topology transition. Once the hot bubbles connect, any new SN 
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will have a much larger V 4 D, and the overlapping will inevitably drive the ISM toward “thermal 
runaway”. Therefore, a steady state with /y,hot ~ 0.6 ±0.1 is better termed as “metastable” due 
to its fragility. 


4-3.3. fy : hot os a function of (n, S) 


Fig. 11 shows how /y hot of the ISM changes when we systematically vary n and S. Each 


simulation is represented by one circle on the h — S diagram, and the colors indicate /y,hot • Each 
data point is obtained by averaging over 20 snapshots evenly sampled during the last 20% of the 
simulation time. From the figure, a lower h and higher S yield a higher /y^ot- As mentioned 
above, those with /y^ot ~ 0.6 ±0.1 lie in the “transitional region” (grey-shaded), where the ISM 
is metastable and can change into “thermal runaway” with any new SN. Below the “transitional 
stripe”, SN do not play a major role in shaping the ISM, and can be significantly under-pressured; 
above that, the hot gas has a “connecting” topology and is in “thermal runaway” state, with hot 
gas dominating the volume and an increasing pressure. 


A simplest estimate of the /y hot is 


47r. 


/v,hot ~ g (<5^cool) ^Close'S*- 


(27) 


(f>Rcooi is the size of the hot bubble, where cf is about 1.6-1.8 (Section 3.1). f c i OS e is the time for a 
SNR to close up: 

^close ~ fiRcool/O s . (28) 

Therefore, 


r mu ^4f 10km / s v S 

Jv.hot ~ U.ll( —) (- Hitt,-=3™ 

1.7 c s 100 kpc d Myr 


■^SN 


- _)(_ 7 _)-!■<*( ^ ) 

3 Mvr -1 1 cm 3 v 10 51 erg 


\ 1.16 


(29) 


Eq. 27 29 only hold true when /y^ot 1$ 0.6, i.e. SNRs do not overlap and thus thermal runaway 


does not come into operation. According to Eq. 29 given an n, /y.hot should be proportional to S 


This seems to be (very roughly) the case for hl_S'50, hl_5100, hl_5200. However, the simulations 
show that for a higher density medium, the dependence of f\\\ m \ : on S is sublinear. For example, 
when h = 3 cm -3 , a factor of 10 increase of S - from 10 3 to 10 4 kpc -3 Myr -1 - only leads to a 
factor of ~ 2 increase in /y,hot- Another way to look at it is to plot the corresponding curve 
of /y,hot = 0.6 using Eq. 29 (black dotted line on Fig. 11). For n > 1cm -3 , the simulations 


with/V’hot ~ 0.6 ±0.1 (grey-shaded) have an S above the simple theoretical expectation at a given 
h; whereas for n < 1 cm -3 , the critical S to achieve a HDMM has a weaker dependence on n than 
expected. Overall, thermal runaway appears harder to achieve in dense gas. 


We argue that one key factor to explain the discrepancy is the photoelectric heating. Recall 
that we have assumed r pe oc S. If we consider a medium without SN (but with normal cooling and 
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PEH), then the warm and cold phases would co-exist with the same pressure, each in its thermally 
stable state. If the mean density of the ISM is fixed, a higher r pe means more mass in the warm 
phase. Note that the cold medium usually occupies a small volume whereas the warm phase 
permeate the space, so more mass in the warm phase suggests that a SNR will travel effectively 
in a “denser” medium, in which it would have a smaller R coo i- This would have a substantial 
negative effect on /v,hot when we consider multiple SNRs, since /v,hot depends on R coo \ sensitively 
as /v,hot oc R* ool oc n _L68 (Eq. 


27 


and 


28). 


To test this idea, we take the run hl0_510000, and reduce its original PEH rate r pe by a 
factor of 10 and 100, respectively. The latter case corresponds to a PEH rate similar to the solar 
neighbourhood. The comparison is shown in Fig. [12} 


The effect of r pe is obvious. Hot gas occupies more volume and has higher pressure for lower 
PEH rates. For the lowest PEH run, the ISM is clearly in a thermal runaway state, since the hot 
gas occupies most volume and determines the pressure. By ”determining the pressure” we mean 
that, if no SN explodes, the pressure of a thermally stable two-phase medium would be about 
3 x 10 3 cm _3 K for a PEH rate of 1.4 x 10~ 26 erg/s, but now the ISM pressure is 1-2 orders of 
magnitude higher than the equilibrium pressure. In this case the overlapping SNRs dominates the 
thermal state of the ISM. The mass fraction of the cold gas is also much higher in the low PEH 
run, which is likely due to two effects: a smaller background heating rate that allows more gas in 
the cold phase, and compression from the overlapping SNRs. Interestingly, the volume-weighted 
pressure does not have a monotonic dependence on r pe : the intermediate r pe case has the lowest 
pressure among the three. This exhibits a subtle balance between the hot and warm phases: SN 
marginally overlap but are able to get rid of most of the energy; the pressure is slightly enhanced 
from that of a two-phase medium. PEH plays the leading role in determining the ISM properties 
in the highest r pe case, whereas SN wins in the lowest one. Note that in Section 3.1, we have 


shown that PEH rate does not impact the SN evolution in a uniform medium. It affects the ISM 
through regulating the warm/cold ratio in the gas reservoir. We emphasize that the PEH is a 
critical process in the ISM - it plays an important role in determining the portion of various phases 
and the overall pressure (Wolfire et al.]|1995 2003 ). This may also have strong implications for the 


global structure of the disk and star formation (Tasker 2011). 


We give a simple fitting formula for the critical SN rate S'crit as a function of n, above which 
thermal runaway occurs, i.e. /v,hot (h, Scrit) ~ 0.6: 


S cri t = 200( — _ 3 ) k ( ) kpc 3 Myr 1 , 

1 cm J 10 01 erg 


(30) 


where k = (1.2, 2.7) for n < 1, > 1cm 3 , respectively. We have shown this with a grey solid line 
in Fig. [TTJ 
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4-3.4- Pressure as a function of (n,S ) 


Fig. 13 is similar to Fig. 11 


but the color-coding indicates the volume-weighted pressure 
Pv/^B for each simulation. The pattern of P(n,S ) is clear. For a fixed n, pressure increases 
with S, which is expected since higher S means higher heating rate. The scaling is nearly linear. 
The dependence of Py on n is less sensitive, and the correlation is mostly negative, at least for 
S < 10 4 kpc -3 Myr -1 . 

Can we understand Py as a function of (n, S) more quantitatively? We have seen that it 
is not easy to obtain a very accurate /v,hot from a simple analytic argument, especially for the 
thermal runaway case when /y,hot keeps rising. On the other hand, if /v,hot is known, can we 
predict pressure? In fact, if we neglect SN and only consider a two-phase medium, one can predict 
quantities such as temperature, pressure, density of each phase for a thermal equilibrium condition, 


given a cooling curve and a heating term (e.g. Draine 1978) (note that under some circumstances, 


there may only be one phase). If the two phases coexist, they will reach pressure equilibrium. The 
equilibrium solution can be shown by a phase diagram, which is the correlation between any two 


14 


of the three thermodynamic variables, density, temperature and pressure. We show in Fig. 
examples of the equilibrium P eq — n relation, given our cooling curve and a variety of PEH rates. 
The equilibrium solution is thermally stable for the part of the curve with positive slope. The 
impact of different PEH rates on the equilibrium pressure of the neutral phase has been explored 


by Wolfire et al. (1995, 2003). Our simple prescription predicts a very similar P eq — n relation to 


their more detailed model. Take the curve with the MW PEH rate for an instance (the third curve 
from the bottom): the unstable solution corresponds to roughly in between 0.4-15 cm -3 , and most 
part outside is stable (except a second dip between 15 — 30 cm -3 ). When the two phases coexist, the 
warm phase has n Wieq ~ 0.4 cm -3 and T ~ 10 4 K, and the cold has n Cjeq ~ 15 cm -3 and T ~ 250 
K. Another way to look at this relation is that, if a mean gas density n is in the range [n W)eq , n c , eq ], 
the medium will settle into two phases with pressure P eq /kB ~ 3 x 10 3 cm -3 K; if n < n Wjeq , the 
medium will be all warm and P oc n; if n > n Cieq , the gas is all cold with roughly P oc h 0 ' 5 


When we consider SN and the hot phase, the hot component can occupy a significant volume 
but usually contains negligible mass, so the mean density for the warm/cold phase is simply n /(1 — 
/v,hot)- Therefore, given an n, PEH rate and /y .hot , the pressure can be predicted in the following 
way. The equilibrium densities of warm and cold phase scale linearly with PEH: 


n 


0.4 cm 3 ( 


w,e q 


pe 


n, 


15 cm 3 ( 


pe 


' 1.4 x 10 -26 erg/s y ’ c,Lq " v 1.4 x 10 -26 erg/s y 
Then the pressure is (assuming the warm phase has a constant temperature 10 1 K): 


(31) 
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(32) 
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In Table[3]we show the ratio between our analytic calculation and the simulation results Pv,pred/Pv■ 
They agree within a factor of 3 generally. Note that in the thermal runaway case, the hot gas 
pressure is higher than that in the warm and cold, but the ratio is on the order of unity. 

In cosmological simulations where the resolution in space and time is low, the multi-phase 
structure of the ISM and individual SN are often not resolved, and the collective SN feedback may 
be modelled as a sub-grid effective equation of state, if the the cell size is much larger than the 

can 
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spatial domain of a single SN, and the timestep is longer than tb (Eq. 21). Although Eq. 
predict the pressure, it has an extra parameter /v,hot which cannot be easily obtained. We are 
then tempted to fit Py directly as a function of n and S. An advantage is that from Fig. [13], the 


pattern of Py(n,S) seems clear and simple and power-law-like. We thus fit the output data using 
the least squares method, which yields: 


Py/^B = 10 ' 


,3.47=1=0.05_-3 


cm 3 K( 


\0.87=1=0.05 


100 kpc -3 Myr 


,.-i 


^ 1 cm -3 


-0.33=1=0.07 


(33) 


The indices show that Py scales almost linearly with S, while it is somewhat less dependant on 
n. The fitting formula gives a reasonable prediction of pressure for h < 30cm~ 3 . We caution that 
the above fitting formula is time-independent. It is fine for an ISM that reaches equilibrium, but 
for those thermal runaway ones, the pressure is not constant but is still rising at the end of the 
simulation. To generalize the result to the thermal runaway cases, we extend Eq. [33] by adding a 
factor of t/t R im to allow for time evolution. Just to remind readers, f s i m = 2Mb, and tb is given in 


Eq. 21 as a function of (n,S). Thus, 


2.9 x 10 3 cm _3 K( 


\0.87 


0.33 


Py/k B ={ ' 10 ° kpC Myr 

' 9.6 x 10 3 cm- 3 K( s 


(uE=*) 

1.87 ( n \—1.59 


_U.87f n \ 

1000 kpc -3 Myr” 1 ' ' 1 cm- 3 > 


if S < S C nt, 

* ) if Sent- 


(34) 


lOMyr 


S cr it is given in Eq. [30] which only depends on n and Esn- We emphasize that the above model 
holds only under the following conditions: (1) the spatial domain is at least a few times the size 
of a single SNR, and (2) the elapsed time t is on the order of a few tb- If t <C tb, the ISM is not 
fully impacted by SN, and if t 3> tb and the ISM is in the thermal runaway state, the gas will 
probably expand due to the increasing pressure. In cosmological simulations, for a “star-forming” 
cell, given its gas density, S, and a time-interval, we can use Eq. [34] to model the SN-modihed 
pressure, provided that the above conditions are satisfied. Note that we do not take into account 
the “turbulence pressure” here. As we will show in the next subsection, the Mach number is close 
to unity, which means that the pressure from the random motions is similar to that of the thermal 
pressure (Joung et al.| 21)1)9). In practice, a factor of 2 can be added to Eq. 34 to make it the total 
pressure. Finally, an interesting coincidence is that the box size and duration of the simulations 
are roughly on the same orders of the cell size (~ 100 pc) and timestep (~ 10 Myr) of cosmological 
simulations, respectively. 
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4-3.5. Other output parameters 


Table [3] lists the input/output parameters for all simulations in this section. For each ISM 
phase, we show the mass fraction f m , volume fraction fy and volume-weighted Mach number My. 
Again, the numbers are obtained by averaging over 20 snapshots evenly sampled during the last 
20% of the simulation time. The Mach number is defined as 


M = v/c s: (35) 

where v and c s are the local velocity (in the simulation box frame) and the local sound speed, 
respectively, and 

c s = y/^ksT/p, (36) 

where the adiabatic index 7 = 5/3, and the molecular weight p = 1.22, 0.588 mjj for T < 10 4 , 
> 10 4 K gas, respectively. 


The order of unity value of M may not be surprising, since the blast waves induced by SN 
explosions would dissipate until they decay to sound waves. Intriguingly, however, the deviation of 
M. is fairly small for the three phases, despite the large parameter space (n, S) we have explored. 
We find My = 0.5 ± 0.2, 1.2 ± 0.3, 2.3 ± 0.9 for the hot, warm and cold phase, respectively. 
The corresponding mass-weighted A4 m (not listed in the table) are 0.5 ± 0.3, 1.0 ± 0.2, 2.6 ± 1.2, 
respectively. The warmer phases have consistently lower M. The Mach numbers of the cold phase 
lie within the observational range M ~ 1—4, as derived from the thermal pressure distribution in the 
solar neighbourhood (Jenkins & Tripp 2011). However, we should caution that the cold clouds are 
usually resolved by < 10 grids in radius, and the velocity dispersion of the gas might be suppressed 
due to the limited numerical resolution. Finally, for a given density, My for the warm/cold phase 
seems to correlate (very) modestly with the S for either h ^ 10 cm ” 3 or h ^ 0.1 cm” 3 . Note that for 
the former n range, the ISM is in a steady state (or metastable), and for the latter, thermal runaway 
(or metastable). Observationally, the velocity dispersion of Ha may have a mild correlation with 
star formation rate (e.g. Green et al.||2014 Arribas et al.|[2014). Whether the velocity dispersion of 


HI has such a correlation is not clear; evidence for either side has appeared in literature (Tamburro 


et al. 2009 Rogers Sz Pittard 2013, Stilp et al. 2013), but the correlation is modest at most, if 


any. Our results are broadly consistent with the observations. However, here we are considering 
only SN as the turbulence driver; other processes, such as gravitational and/or magneto-rotational 
instability, can also play roles in driving the turbulence. Also note that for runs with larger h, 
we use higher resolutions and smaller boxes. Thus the finer structures are increasingly resolved, 
especially for the dense phase. For the lower n runs, the resolution of the inner structure of the 
cold clumps may be limited by the finite cell size. 


In Table [3] we have also shown the fraction of the kinetic energy out of the total energy in 
the simulation domain, fy k- This is obtained by averaging over 30 snapshots for the last 3% of 
the simulation time. fy,k is in the range of 15% - 50%, and the value is smaller for lower n. This 
is consistent with that the denser media cool more efficiently, thus leaving a larger fraction of the 
energy in the kinetic form. 
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We define a “heating time scale” as 


_ e tot,V _ e tot,V 

esN S ■ Esn 


(37) 


where etot is the volume-weighted energy density of the ISM at the end of the simulation, esN is the 
time-averaged SN heating rate, and FgN = 10 51 erg. If the ISM reaches the thermal equilibrium, 
egN is equal to the time-averaged cooling rate e coo i (the other heating source, PEH, is an order 
of magnitude weaker than the SN overall), so fheat is approximately the cooling timescale. If the 
ISM is in the thermal runaway state, fheatAsim ~ 1 — e coo i/esN, where t s j m is the duration of the 
simulation. From the table, fheat is generally smaller than 1 Myr, much smaller than £ s ; m , indicating 
that cooling is efficient in the ISM. Even in a thermal runaway state, e coo i/esN is very close to unity. 
For extreme cases, such as h0035_5100, fheat/4im ~ 0.2, which implies that on average 20% of the 
energy from each SN is retained in the ISM, instead of being radiated away. For the majority of 
the thermal runaway media, the retained fraction of energy per SN is about 1-5%. 


Discussion 


5.1. Randomness of SN position on simulation box scale 


In the simulations described in Section 4, we have kept the positions of SN random. In this 
section we discuss how good this assumption is. Type la SN, which are from old stellar populations, 
are not correlated with the gas properties, and are thus close to random. Core collapse SN are more 
complicated. They descend from OB stars, and are mostly associated with star-forming clusters, 
which have abundant dense molecular clouds. Recent work have shown that the SN position relative 
to the dense clumps has a significant impact on the feedback. Iffrig & Hennebelle (2015) found 


that a SN is most disruptive to a molecular cloud when it explodes inside; a SN outside has only 
moderate effect on the cloud, but the feedback influence is felt much further in the low-density 


inter-cloud medium. Considering multiple SN, Gatto et al. (2015) discovered in a periodic box 


study that a random positioning leads to thermal runaway more easily than when the SN mostly 


explode within the clumps. Hennebelle Sz Iffrig (2014) used a set-up with stratified medium and 


found that random SN lead to a higher gas density and pressure at high galactic latitude - closer to 
what has been observed. The random positioning leaves most SN to explode in a rarefied medium, 
which enables an effective feedback on larger scales. 

To assess the SN position relative to the dense clumps, we emphasize that the non-negligible 
velocities of the OB stars should be taken into account. We note that some other effects can also 
lead core collapse SN to explode in a low-density environment. For example, molecular clouds 
disperse on a timescale of about 20 Myr, shorter than the life time of some SN progenitors, so even 
those in situ may explode in a rarefied gas; recent works (e.g. Bressert et al.||2012 |Oey et al.| 2013) 
suggested that some OB stars are indeed born in the field, rather than in crowded clusters. But 
here we simply focus on the effect of the OB velocities. In particular, it would be very interesting 
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to know how far OB stars can migrate before they explode as SN. To know the answer helps other 
important questions: (1) What fraction of core collapse SN explode inside molecular clouds versus 
outside? (2) Is the random positioning of SN assumed in our simulations a good approximation? 

To calculate the displacement distribution of core collapse SN progenitors, we do the following 
simple Monte Carlo simulation: We adopt the distribution of space velocities of OB stars compiled 


by Stone (1991), which is a sum of two Maxwellians with a = 7.7 km/s and 28.2 krn/s, respectively, 
and 46% of O stars and 10% of B stars belong to the high-velocity group. More recent work by de 
Wit et al. (2005) confirmed that roughly half of the O stars are “runaways”. Each core collapse 


SN progenitor is randomly assigned a mass in the range 8M 0 < M* < 50 Mq according to the 
initial mass function dN/dM * oc M~ 2,3 (Kroupa 2001, Chabrier 2003), a speed according to the 


distribution mentioned above (if M* > 16 Mq, the progenitor is classified as an O star, otherwise a 
B star), and a random velocity direction. The life time of the progenitor is related to its mass by 
tiife = 10 10 year(M*/ Mq) -2 " 5 . 


The displacement distribution is shown in Fig. 15, In 3D space (left panel), about 96%, 82%, 
70% of SN progenitors migrate away from their birth sites by distances greater than 10, 50,100 pc, 
respectively. The giant molecular clouds appear filamentary (e.g. Molinari et al. 2010), and the size 
is on the order of 100 pc in the long dimension and an order of magnitude smaller on the shorter 
ones. Our result suggests that roughly < 5% of SN would explode within these clouds. Those SN 
can cause the clouds to disperse and thus suppress star formation. On larger scales, the external 
gravitational field from the dark matter halo and stellar disk may preferentially confine the motion 
of the progenitor stars to a certain direction, so it makes more sense to calculate the projected 


ID displacement distribution (right panel of Fig. 15). Approximately 17%, 8.4%, 2.5% of the SN 


progenitors have displacements greater than 0.5, 1, 2 kpc, respectively. A significant fraction of 
runaway OBs are therefore likely to migrate out of the star formation regions, and explode in the 
low-density part of a galaxy, such as halo and inter-spiral arm regions, and may drive a wind there. 

The sizes of our simulation domain range from a few tens to a few hundred parsec, and we want 
to know how the runaway stars affect the randomness of the SN position on this scale. We conduct 
the following comparison experiments: we run two additional simulations to hl_S'200 (“Random”) 
- (1) the probability of where a SN progenitor is born is proportional to local gas density p 1 ' 5 , 
and each SN is randomly assigned a displacement according to the displacement distribution in 
the left panel of Fig. 15 (“HD+runaway”, where HD stands for high-density); (2) same as (1) but 
no displacement, that is, SN explode where their progenitors are born (“HD”). All other input 
parameters, such as n, S, the box size and resolution, remain the same. 


The results are shown in Fig. 16 the properties of the ISM produced by “HD+runaway” fall in- 
between those in the other two cases. In particular, displaced SN yield very similar volume/mass 
of each gas phase to random SN. The “HD” run produce more warm gas in mass and a much 
smaller volume fraction of the hot medium. This is understandable since the density peaks in the 
ISM consist mostly of the cold medium, and the energy released by a SN converts the cold medium 
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to warm. The strong cooling in those dense clouds restrains the power of SN (i? C ooi oc n” 0 ' 42 ), so 
they cannot impact the ISM on larger scales. This explains the smaller hot gas fraction and lower 
pressure. Those results are qualitatively consistent with Gatto et al. (2015). We conclude that 


based on the observed OB star velocities, the positions of core collapse SN are nearly random on 
scales < 150 pc. 


5.2. Connecting simulation data to real galaxy environments 


We have plotted ( n,S ) encountered in the MW disk to halo on Fig. 11 (blue solid line), as 
described in Section [4.2| For the disk average hl_S'200, the ISM is in the “transitional stripe”. Of 


course, the disk average does not take into account the clustering effect of the SN. OB stars are 
born mainly in OB associations, and eventually 60% of core collapse SN appear “grouped”. Since 
the disk average S makes the ISM in the transitional region, it is not surprising that any grouping 
in space and/or time can lead to thermal runaway, assuming the same gas mean density. Indeed, 
if SN explode within a sphere of R ~ 100 pc every 0.3 Myr as for a super-bubble (Mac Low & 
McCray 1988), which corresponds to S ~ 800 kpc~ 3 Myr -1 , it is well within the thermal runaway 
regime for h = 1cm -3 (Fig. 11). 


As one goes up to the halo, the mean density is smaller, and the SN rate is such that /v,hot is 
at first larger and the ISM goes into the thermal runaway region. This suggests that at least part 
of the halo is not in hydrostatic equilibrium. Above a certain height, fy hot drops again, which may 
suggest that the gaseous halo is convective, but a conclusion cannot be made without a simulation 


with stratified medium and an external gravitational held (Joung &; Mac Low 2006; Gent et al. 
2013 Creasey et al.||2013 ). 


We have also plotted the scaling relation of S and h as suggested by the Schmidt star formation 
relation S oc n 1-5 (green dashed line). The relation is normalized such that the line passes through 
the MW average hl_S’200. For data points that lie around this relation, /v,hot is smaller for higher 
densities. For n = 30 cm” 3 , the ISM around the Schmidt relation has a /y,hot °f on ly a few percent. 
Surprisingly, for ISM in star-burst regimes, little hot gas exists and the medium is thermally stable. 
Observationally, however, powerful winds are ubiquitously associated with such high star formation 


rates, and the mass loading (= M out fl ow /MsF) is on the order of unity or even higher (e.g. Heckman 
et al.|2000 Steidel et al.|2004 ). Note that we have found the PEH rate to be a very important factor 
in determining the thermal state of the ISM. If the PEH rates we have adopted are reasonable, then 
the data seem to suggest that one needs other mechanism(s) to drive a wind for the high density 
regions. For example, “clustering” of SN may work, although a factor of 10 or more increase in S 
is needed for n > 3 cm” 3 above the average value to have a thermally-driven wind. Another factor 
is the pre-SN feedback, such as photo-ionization, radiation pressure and winds from massive stars, 


which can create low density tunnels, facilitating SN energy to leak out (e.g. Rogers & Pittard 


2013). Alternatively, runaway OB stars, which can easily migrate more than a hundred parsecs, 


may lead to a significant fraction of core collapse SN exploding outside the immediate high density 
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SF clouds. For a low density medium n < lcm' 3 , the critical SN rate for the ISM to have a 
thermal runaway is much less stringent. Another possible mechanism for wind launching is by the 
cosmic rays. As the cosmic rays diffuse out, the pressure gradient exerts accelerating force on the 
baryonic gas. Recent simulations have shown that this mechanism is promising to drive winds with 
a reasonable mass loading (Uhlig et al.||2012 Hanasz et al.|[20l3 Salem & Bryan 2014). 


5.3. Periodic box 


In the experiments where multiple SN shape the ISM, we have applied periodic boundary 
condition and initially uniform medium. We briefly discuss their applicability here. 

In a periodic box, the total mass is fixed, and the gas is not allowed to expand beyond the 
simulation domain. This means that such a box is good at capturing processes that have scales 
smaller than its size, which do not result in net ingoing/outgoing flux on the outer boundary. For 
our experiments, this is the case for those reaching steady states. For the thermal runaway cases, 
the pressure keeps rising and the box is accumulating energy. In real galactic environment, the hot 
gas is subject to expansion if the ISM pressure around is not sufficient to stop it. The expansion 
can lower the pressure and /y,hot (de Avillez & Breitschwerdt 2004). Such an evolution thus cannot 
be faithfully captured by a periodic box. 


On the other hand, we want to emphasize that the onset criterion for the thermal runaway, 
i.e. when the SNRs overlap, as given in S’ c irt(h), should not depend on whether the box is open or 
periodic - it only depends on the local density of the gas. Furthermore, even initially the outgoing 
flow is in an explosive state, it may adjust so that each patch of the flow comes to an equilibrium 
state. The length scale of the galactic winds is ~1-10 kpc, much larger than R C ooi of a single SN, 
and the dynamical time scale is ~ Gyr, much longer than what takes to reach a steady state. So 
when the gas on large scale is in a steady state, the effect of SN may be regarded as “microphysics” 
and incorporated in the effective equation of state in a cosmological simulation. 


A related issue is the density distribution in the box. In our simulations we keep gas uniform 
as the initial condition. We note that in real galaxies, especially disk galaxies, the ISM is stratified 
in the vertical direction. Therefore, the uniform assumption only holds when the processes we 
study have scales smaller than the scale height of the gas. The size of our simulation domain, 
which is 6i? coo i) is generally smaller than the gas scale height in the MW, given a mean ISM 
density. For example, the thickness of the neutral gas layer is about 180 pc with a mean density 
~ 1 cm -3 , the warm ionized medium has a scale height of about 1830 pc with the density 0.014 cm -3 
(Gaensler et al. 2008). For comparison, the box size is 142.2pc for n = lcm' 3 , and 585.6pc for 
h = 0.035 cm -3 . So the non-stratified medium is probably not a bad assumption. That said, the 
gas scale height in extragalactic systems are usually not well-constrained, and SN themselves may 
play an important role in determining it (e.g. Shetty &; Ostriker 2012). 


In disk galaxies, clustered SNR in a thermal runaway state can form “super-bubbles” that 
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may break out of the disk and transfer materials into the halo (Tomisaka & Ikeuchi 1986; Mac 


Low & McCray|[l988;). The outflowing gas forms galactic winds and/or fountains that regulate the 


circum-galactic environment (|Heckman et al.|2000 Steidel et al.|2004). These can only be captured 


in a larger box with more realistic set-up. For example, the medium is stratified, which has outflow 


outer boundary conditions in the vertical direction and external gravitational field (e.g. Joung & 


Mac Low 2006 Creasey et al. 2013). We will postpone the study of SN feedback in such cases to 


future work. 


5.4. Thermal conduction and magnetic fields 


Thermal conduction can be an important process in a multiphase medium to drive mass and 


energy transfer among different phases. In the model of McKee & Ostriker (1977), evaporation 


through thermal conduction plays a major role in producing the warm phase from the dense, 
cold clumps in the post-shock region. Very small clouds bathed in a hot medium can be quickly 
evaporated. As mentioned earlier, thermal conduction is not included in our simulation. The warm 
phase is maintained mainly through the balance of the shock heating/PEH and cooling. We did 
some calculations to estimate the role of conduction for our simulations. For the multiphase ISM 


models we adopt in Section 3.2, conduction actually causes hot gas to condense on the clouds (Cowie 


& McKee 1977), and the evaporation is only important when very close to the explosion center, 


where the temperature is high and clouds are shock-stripped to small cold cores. The difference 


between our simulation and McKee & Ostriker (2007) arises from that our clouds, which are created 


and destroyed by SN interactions, are larger by a factor of a few than what was adopted in the 


McKee &; Ostriker (2007) model, which was inspired by the HI observation at that time. Smaller 


clouds would result in a more effective mixing between the different phases in a number of ways, 
e.g. Kelvin-Hemholtz instability, and photo-evaporation by massive stars (McKee et al. 1984). To 


fully evaluate the role of conduction, one also has to consider the magnetic field, which can cause 
significant anisotropy. Conduction is practically inhibited perpendicular to the field lines. When 


a blast wave travels in a HDMM, it tends to wrap the field lines around the clouds (Semenov & 


Bernikov 1980; Dursi & Pfrommer 2008), which suggests that the conduction between clouds and 


the hot medium may not be very important. The actual impact of conduction considering magnetic 
field is not clear and is definitely worth further study. 

The magnetic field itself can also be important dynamically. Local observations show that the 
magnetic pressure is a few times larger than gas thermal pressure (Heiles & Crutcher 2005). Slavin 


Sz Cox (1992) has found that in a uniform medium, the magnetic pressure at late stages of the SNR 


can thicken the shell and thus make the hot bubble smaller. This may change the critical SN rate 
for thermal runway at a given density. A MHD simulation is needed to quantify the role of the 


magnetic field in shaping the ISM (e.g. Shin et al. 2008 Hill et al. 2012). 






























































6. Conclusion 


In this paper we explore the interplay between SN and the three-phase ISM. We conduct three 
sets of experiments to compare the evolution of a SNR in a uniform and a hot-dominant multiphase 
medium, and to allow SN to self-consistently shape the ISM under various environments. The main 
results are summarized as follows: 


1. For a uniform medium, the evolution of a SNR has four well-defined stages after the initial 
free expansion: Sedov-Talor, (transient) “pressure-driven snowplow”, momentum-conserving and 
close-up. Thermal energy drops precipitously once the shock becomes radiative as E t h oc R“ 6 ~ -8 , 
while kinetic energy drops later and more slowly as E^ oc R p 3 . The R p — t relation and the 
energy/momentum feedback vs R p can be well described by piece-wise power-law functions. The 
power-law indices vary little over a wide range of gas densities and photoelectric heating rates 
(Table 1). 


2. In a hot-dominated multiphase medium, a SNR evolves very differently from the uniform 
case or a two-phase medium. There are no distinct stages of evolution. The blast wave initially 
interacts with the clouds, but later travels much faster and mainly in between the clouds (Eq. 13). 
The SNR thus sweeps up much less mass, therefore developing less momentum at a given R p (Fig. 
[6]). A radiative-efficient thin shell never forms at the blast wave front, and the overall cooling rate 
is suppressed. The hot phase helps to retain the SN energy and facilitates energy transfer to much 
larger scales. The spatial-temporal domain a SN is enlarged by a factor of > 10 2 ' 5 in a HDMM. 


3. We study the ISM self-consistently shaped by multiple SN, covering a wide range of SN 
frequency S and gas average density h. We find that the classic view of the three-phase ISM shaped 
by SN, as envisioned by McKee Sz Ostriker (1977), needs to be extended, depending on n and S. 
The main differences we have found are the following: First, a certain phase can be absent, or exist 
only very briefly. Second, equilibrium solutions may not exist when SN explode at a sufficiently 
high rate, and the ISM undergoes thermal runaway. Third, the pressure of different phases may 
not be in equilibrium; the hot phase can be over- or under-pressured. 


4. We find that /v,hot ~ 0.6 ±0.1 marks the transition of ISM from a steady state to thermal 
runaway, along with the change of hot gas topology. When /y^ot 0.6 ± 0.1, SNRs evolve as 
individual bubbles; if /y^ot can reach 0.6 ±0.1, the hot bubbles connect to each other. The 
connecting topology makes any new SNR have a larger spatial-temporal domain of influence, which 
reinforces the overlapping of the blast waves. As a result, the overall heating dominates over cooling 
and the ISM undergoes thermal runaway. An ISM with fvhot ~ 0.5 — 0.7 is in a metastable state, 
and may transition into thermal runaway. 


5. The PEH has a surprisingly strong impact on /y^ot- For a fixed (h,5), /v,hot decreases 
as the PEH rate increases. In particular, for n > 3 cm -3 , assuming that the PEH is proportional 
to S, f\r t hot is confined to < 0.6 for a wide range of S, so that the medium is stable and does 
not form a wind (Fig. 11). The critical SN rates for the onset of the thermal runaway is roughly 
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5'crit = 200 (n/1 cm 3 ) fc (-E’sN/10 51 erg) 1 kpc 3 Myr 3 , where k = (1.2, 2.7) for n < 1 and > 1 cm 3 , 
respectively. 


6. The pressure of the ISM relates to /y^ot; mean gas density n and photoelectric heating 


rate T pe through Eq. 31 32 The results do not depend on whether the ISM is in a steady state or 
thermal runaway. A simple fitting formula Pv{n, S, t) can be used as an effective equation of state 
considering SN feedback: 


Pv/^b = 


2.9 x 


103cm 3K Uok P A ly ,-± 87 <±±> 033 if s< Sci., 

9.6 X ltf if 


7. The local Mach numbers A4 of the three ISM phases show surprisingly small variations, 
despite the 5 orders of magnitude span of (h, 5). We find that My ~ 0.5 ±0.2, 1.2 ±0.3, 2.3 ±0.9 
for the hot, warm and cold phase, respectively. 

8. We calculate the displacement distribution of the core collapse SN from the observed 
velocities of OB stars, which shows that on scales < 150 pc, SN explode almost at random positions, 
and nearly 10% of OB stars can migrate > 1 kpc (Fig. [15] ). The latter runaway OBs have a chance 
to migrate into low-density regions where it is easier to drive a wind. 
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Fig. 1.— Evolution of a SNR of initial energy 10 51 erg in a uniform medium with n = lcm~ 3 . 
Snapshots are taken at t = 1.5 X 10 4 , 4.6 X 10 4 , 9 X 10 4 , 2.9 X 10 5 , 1.1 x 10 6 , 2.6 X 10 6 , 4.2 x 10 6 years, 
respectively. Numerical resolution is 0.75 pc. Vertical axes show spherically-averaged quantities, 
and horizontal axes are radii. Dashed lines in the upper-left panel show the positions of the 
shock as predicted by the ST solution at the above timesteps. Pressure of the inner bubble drops 
precipitously once the thin shell forms, to even below that of the ISM. At later stage, the shell 
becomes thicker and keeps moving outward to R > 80 pc, while the the hot bubble (T > 2 x 10 5 
K) barely reaches 42 pc in radius. Velocity is shown in absolute value. 
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Fig. 3.— Slices of density, temperature, pressure and velocity of a 10 51 erg SNR in ISM model “fvhO.GOa” with n = 1 cm 3 (see 
Table [ 2 ] for full parameters). Snapshots are taken at 1.9 x 10 4 year (upper panels), 1.1 x 10° year (middle) and 2.6 x 10 5 year 
(lower), respectively. The velocity plot has a floor of 0.1 km/s for log-scale display. Note that the first two rows have zoomed-in 
views. 
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Fig. 5.— Snapshots for a 10 51 erg SNR evolving in ISM models “fvh0.66ss” (upper row), “fvh0.90a” (middle), “fvh0.82ss” (lower), 
respectively, at 3 x 10 4 year. See Section 3.2.1 for model details. The velocity plot has a floor of 0.1 km/s for log-scale display. 
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Fig. 6.— SN feedback in the 4 multiphase ISM models described in Section 3.2 vs a uniform 
medium with the same mean density n = lcm -3 . Line styles and the corresponding models are 
given in the legend in panel (c). (a) R p vs time; (b) momentum vs R p \ (c) E t h vs time; (d) E t h vs 
R p , (e) Ek vs time; (f) E^ vs R p . In (b), only results of “fvh0.60a” and “fvh0.82a” are shown, and 
the thin and thick lines denote the total momentum and the radial component, respectively. The 
black dotted lines in each panel indicate the fitting formula presented in Section 3.2.3| 



































Fig. 7.— Slices of density, temperature, pressure and velocity of the ISM from simulations with varying mean density n and 
SN rate S. Upper row: hl_£200 (solar neighbourhood) at t = 104 Myr; middle: h0.3_5100 at t = 46 Myr; lower: hl0_S3000 at 
t = 135 Myr. Note that the box sizes vary. See Section [4. l| for model details. Most volume is shared by the three typical phases 
with T ~ 10 6 ,10 4 ,10 2 K, respectively. Depending on n and 5, a certain phase can be missing. The clumping of gas varies with 
n and S. 
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Fig. 8.— Evolution of the three ISM phases for nl_S'200. Up to down, left to right: volume- 
weighted pressure Py/ks , Mach number, volume fraction, mass fraction. Red triangles, yellow 
diamonds and blue circles denote the cold (T < 10 3 K), warm (10 3 K < T < 2 x 10° K ) and hot 
(T > 2 x 10 5 K) phases, respectively. The medium reaches a steady state and the three phases are 
in rough pressure equilibrium. 
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Fig. 9.— Volume-weighted pressure Py of the three phases versus time for runs h0.1_S'50 (left) 
and hl0_S'3000 (right). For n0.1_S50, the hot gas is over-pressured and the ISM is in a “thermal 
runaway” state. For nl0_S'3000, the hot phase is only present briefly after SN explosions and is 
significantly under-pressured. 



Time/Myr Time/Myr 

Fig. 10.— Transition from a metastable state to “thermal runaway” for nl_S'200. Initially, the ISM 
reaches the apparent steady state, when roughly half the volume is hot. The transition happens at 
t ~ 55 Myr, when the hot gas quickly dominates pressure and > 90% of the volume, and pressure 
of all phases rises progressively. 
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Fig. 11.— Volume fraction of the hot (T > 2 x 10 5 K) gas /y^ot as a function of the mean gas 
density n and SN rate S. Each circle represents one simulation. Colors show the /y^ot • The 


blue solid line indicates S as a function of h from MW mid-plane to halo (see Section 4.2 for 
a description). The green dashed line indicates S oc n 1 " 5 . (n, S) = (1,200) is roughly for the 
MW-average. The shaded area shows the transitional region, i.e. /v,hot ~ 0.6 ± 0.1 - above it 
the hot gas dominates the volume and the ISM undergoes thermal runaway, and below it the SN 
play a subordinate role and most volume is warm/cold; within the transitional region the ISM is 
metastable. The gray solid line indicates the nominal line to roughly represent above which 


thermal runaway happens (Eq. 27). The black dotted line shows /y,hot = 0.6 from the simple 
analytical expectation Eq. 
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Fig. 13.— 


Same as Fig. 11 


but colors show the volume-weighted pressure Py /fee • 


























- 55 - 



Fig. 14.— Equilibrium pressure as a function of density n for T < 10 1 K gas for a variety of 
photoelectric heating rates. Red to green: 0.1, 0.5, 1, 2, 4, 10, 30, 100, 200 x 1.4E-26 erg/s per H 
atom. The calculation is based on the cooling curve adopted in our simulation. 































- 56 - 



0 < p)d 


o 



(* < p)d 


C/T 

Q 


o 

T—i 


4-2 

Td 


3 

CD 

bjO 

CD 

-4-2 

CD 

OD 


O 




C 


a 

&H 


?—i 

ft 


dJ 

CD 


ft 

ft 


4-2 

+2 





O 

-1-2 





CD 

bO 


8 

s 


a 


d 

ft! 

CD 

Oh 

o 

4-2 



f—1 

o 

a 

ft 

o 

Th 

CD 

T-1 

CD 

ft 

R 

CD 

4-2 

cd 

"D 

O 

ft 


4-2 


CD 

0) 

CD 

r_ D 

_> 

o 

o 

4-2 


p 

cd 

R 




CD 

CD 

T3 

ft 


CD 

4-2 

£ 

Id 

f—1 

o 

CO 

buO 


CD 


CN 

l/J 

a 

a 


p^h 

|> 


o 

a 

cd 

.2 

’-4-2 



CD 

CD 

S-H 

o 

£ 

ft 

0) 

CO 

CD 

<4-1 

CD 

CD 

CD 

o 

a 

CO 

4-2 

cd 


3 

' — 

6 

CD 

O 

a 

a 

CD 

ft 

CD 

CD 

T—1 

£ 

S-4 

o 

A? 

d 

CO 

33 

Q 

<4-1 

o 

££ 

00 

CD 

4-2 

cd 

Jh 

.h0 

CO 

CD 

a 

<4-1 



o 

P3 

pp 

rD 

cd 

o 

bn 


4-2 

R 

C/3 

R 

ft 

o 

PS 

cd 

4-2 

m 

2 


-1-2 

CO 

"3 

CO 

cd 

PQ 

O 

.Sl¬ 

id 

CO 

PQ 

o 

<4-4 

O 

6S 

o 


O 

T— 1 

s 

CO 

D 


r 

S-4 

u 

'-4-2 

cd 

CD 

<s 

CD 

O 

£ 

CD 

O 

> 

423 

a 

■ 

0) 

CD 


ft 

a 


-4-2 

CD 

lO 


c> 

T—1 

-4-2 

dS 

bb 

CD 

3 

a 

m 

d 


X5 















Mass fraction Volume fraction . cnn Pressure / kg(cm 3 K) 

-■-' 1 i '-■-■ 1 45001—■-'-- 



o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

l o 

o 

LO 

o 

LO 

o 

LO 


co 

co 

<N 

<N 

T—1 

T—1 




Lime 

SN 

D 

-4-3 

cS 

J 

ctf 

5 

3 -** 

o 

-4-3 


cb 

u 

on, \ 

“HD 

D 

2 

o3 

Q 

ffi 

+ 

’-4-3 

- 




o 

a3 


ex 

O 

a 


bJO ^ 
§ 

n *-■ 


CD 


o 

a 


<D 

N 


O 

o 

o 

CN 


►-O 

-1-3 


£ 

D 

3 

bO 


-1-3 

c5 

D 

03 

O 

3h 

X 

D 

£ 

CO 




D 

- 4-3 


D 

D 

a3 


rt 

a3 


ex 

o 

a3 

r-3 

a 


03 

a3 

D 

£ 

r^H 

CO 


CO 



CO 


n 

a3 

+3) 

3 

co 

D 

&—i 

D 

-4-2 

n 

o 

a 

.2 

’-+3 

cb 

D 

JO 

!z; 

CO 


D 

D 

?—i 
- 4-3 

D 

- 4-3 


I 

ri 

3 

+ 

Q 

ffi 


D 

J-i 

0 

C/2 

X 

D 

Sh 

a 

T3 

D 

- 4-3 

r-a 

bJO 


- 4-3 

’tf2 

Ch 

D 

3 

3 

bJO 


a 

*§ ^ 

H 

c5 


3 ^ 


3 


T—I 

.Sf 

3 


3 

.2 

’-4-3 

D 

a3 

,S-4 


D 

03 

O 

3h 

x 

D 


velocities of their progenitors. The ISM from the “Random” and “HD+runaway” are similar, whereas that from 
much smaller /v,hot and more warm gas. 



































